      MODULE GWFSWIMODULE
        TYPE TSWIDE4
          INTEGER  :: MXITER,NODES,NHALFU,NHALFL,NBWGRD
          INTEGER  :: MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR
          INTEGER  :: NITERDE4,IFREQ,ID4DIM
          INTEGER  :: NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
          REAL     :: ACCLDE4,HCLOSEDE4,DELTL
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   :: IUPPNT
          INTEGER,          ALLOCATABLE, DIMENSION(:,:,:) :: IEQPNT
          REAL,             ALLOCATABLE, DIMENSION(:,:)   :: AU
          REAL,             ALLOCATABLE, DIMENSION(:,:)   :: AL
          REAL,             ALLOCATABLE, DIMENSION(:)     :: D4B
          REAL,             ALLOCATABLE, DIMENSION(:)     :: HDCGDE4
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   :: LRCHDE4
        END TYPE TSWIDE4
        TYPE TSWIPCG
          INTEGER :: MXITER,NODES
          INTEGER :: ITER1,NPCOND,NBPOL,NITER
          REAL    :: ZCLOSEPCG,RCLOSEPCG,RELAXPCG,DAMPPCG
          REAL    :: DAMPPCGT
          INTEGER :: IHCOFADD
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: VPCG
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: SS
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: P
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: HPCG
          REAL,             ALLOCATABLE, DIMENSION(:,:,:) :: CD
          REAL,             ALLOCATABLE, DIMENSION(:,:,:) :: HCSV
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   :: LHCH
          REAL,             ALLOCATABLE, DIMENSION(:)     :: HCHG
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   :: LRCHPCG
          REAL,             ALLOCATABLE, DIMENSION(:)     :: RCHG
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IT1
        END TYPE TSWIPCG
        TYPE TSWIPCGU
          INTEGER  :: MXITER,ITER1C,NPC,NORD,NITERC,NNZC,NIAC
          INTEGER  :: NIAPC,NIWC,NLAPC,NUAPC
          INTEGER  :: NNZAPROC,NIAPROC
          INTEGER  :: NODES
          REAL     :: ZCLOSEPCGU,RCLOSEPCGU
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   :: LHCHPCGU
          INTEGER,          ALLOCATABLE, DIMENSION(:,:)   :: LRCHPCGU
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IT1PCGU
          REAL,             ALLOCATABLE, DIMENSION(:)     :: HCHGPCGU
          REAL,             ALLOCATABLE, DIMENSION(:)     :: RCHGPCGU
          INTEGER,          ALLOCATABLE, DIMENSION(:,:,:) :: NODEC
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)     :: BC
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)     :: XC
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)     :: AC
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)     :: TAPC
          DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)     :: APC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IAC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: JAC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IUC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IXMAP
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IAPC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: JAPC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IRAPC
          INTEGER,          ALLOCATABLE, DIMENSION(:)     :: IAPCMAP
C           WORKING ARRAYS
          INTEGER,         ALLOCATABLE, DIMENSION(:)      :: IWC
          DOUBLEPRECISION, ALLOCATABLE, DIMENSION(:)      :: DC
          DOUBLEPRECISION, ALLOCATABLE, DIMENSION(:)      :: PC
          DOUBLEPRECISION, ALLOCATABLE, DIMENSION(:)      :: QC
          DOUBLEPRECISION, ALLOCATABLE, DIMENSION(:)      :: ZC
          DOUBLEPRECISION, ALLOCATABLE, DIMENSION(:)      :: SCL
        END TYPE TSWIPCGU

        TYPE TSWIOBS
          CHARACTER (LEN=12) :: OBSNAM
          INTEGER :: KLAY
          INTEGER :: IROW
          INTEGER :: JCOL
        END TYPE TSWIOBS

        INTEGER, SAVE, POINTER :: NPLN,ISTRAT,NOPT,NZONES

        INTEGER, SAVE, POINTER :: NADPTFLG
        INTEGER, SAVE, POINTER :: NADPTMX
        INTEGER, SAVE, POINTER :: NADPTMN
        REAL, SAVE, POINTER    :: ADPTFCT
        INTEGER, SAVE, POINTER :: IADPT
        INTEGER, SAVE, POINTER :: IADPTMOD
        REAL, SAVE, POINTER    :: ADPTVAL
        REAL, SAVE, POINTER    :: SWIDELT
        
        INTEGER, SAVE, POINTER :: NOBS
        INTEGER, SAVE, POINTER :: IOBSHEADER

        INTEGER, SAVE, POINTER :: ISWIZT,ISWICB,ISWIOBS
        INTEGER, SAVE, POINTER :: NLAYSWI
C         SOLVER
        INTEGER, SAVE, POINTER :: NSOLVER
!        REAL, SAVE, POINTER    :: ZCLOSE
        INTEGER, SAVE, POINTER :: IPRSOL
        INTEGER, SAVE, POINTER :: MUTSOL
C         SWI PARAMETERS
        REAL, SAVE, POINTER    :: TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA
        INTEGER, SAVE, DIMENSION(:,:), POINTER :: IBO
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWIHCOF
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWISOLVCR
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWISOLVCC
        REAL, SAVE, DIMENSION(:,:), POINTER :: SWISOLVCV
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: ZETA
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: ZETAOLD
!        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: ZETAIT
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: SSZ
        REAL, SAVE, DIMENSION(:), POINTER :: EPS
        REAL, SAVE, DIMENSION(:), POINTER :: NUS
        REAL, SAVE, DIMENSION(:), POINTER :: DELNUS
        REAL, SAVE, DIMENSION(:), POINTER :: NUPLANE
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICR
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICC
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICUMCR
        REAL, SAVE, DIMENSION(:,:,:,:), POINTER :: SWICUMCC
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: NUTOP
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: NUBOT
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: BRHS
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: QZEXTRA
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: QZEXTRACUM
        DOUBLE PRECISION, SAVE, DIMENSION(:,:), POINTER :: DUM
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: RHSFRESH
        REAL, SAVE, DIMENSION(:,:,:), POINTER :: HCOFFRESH
        INTEGER, SAVE, DIMENSION(:,:,:,:), POINTER :: IPLPOS
        INTEGER, SAVE, DIMENSION(:,:,:), POINTER :: IZONENR
        REAL,    SAVE,   POINTER, DIMENSION(:,:,:) ::SC1
C---------STORAGE FOR OBSERVATION DATA
        TYPE (TSWIOBS), SAVE, DIMENSION(:), POINTER :: SWIOBS
C---------STORAGE FOR SOLVERS
        TYPE (TSWIDE4),  SAVE, POINTER :: SWIDE4
        TYPE (TSWIPCG),  SAVE, POINTER :: SWIPCG
        TYPE (TSWIPCGU), SAVE, POINTER :: SWIPCGU

        TYPE GWFSWITYPE
          INTEGER, POINTER :: NPLN,ISTRAT,NOPT,NZONES

          INTEGER, POINTER :: NADPTFLG
          INTEGER, POINTER :: NADPTMX
          INTEGER, POINTER :: NADPTMN
          REAL, POINTER    :: ADPTFCT
          INTEGER, POINTER :: IADPT
          INTEGER, POINTER :: IADPTMOD
          REAL, POINTER    :: ADPTVAL
          REAL, POINTER    :: SWIDELT

          INTEGER, POINTER :: NOBS
          INTEGER, POINTER :: IOBSHEADER

          INTEGER, POINTER :: ISWIZT,ISWICB,ISWIOBS
          INTEGER, POINTER :: NLAYSWI
C           SOLVER
          INTEGER, POINTER :: NSOLVER
!          REAL, POINTER    :: ZCLOSE
          INTEGER, POINTER :: IPRSOL
          INTEGER, POINTER :: MUTSOL
C           SWI PARAMETERS
          REAL, POINTER    :: TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA
          INTEGER, DIMENSION(:,:), POINTER :: IBO
          REAL, DIMENSION(:,:), POINTER :: SWIHCOF
          REAL, DIMENSION(:,:), POINTER :: SWISOLVCR
          REAL, DIMENSION(:,:), POINTER :: SWISOLVCC
          REAL, DIMENSION(:,:), POINTER :: SWISOLVCV
          REAL, DIMENSION(:,:,:,:), POINTER :: ZETA
          REAL, DIMENSION(:,:,:,:), POINTER :: ZETAOLD
!          REAL, DIMENSION(:,:,:,:), POINTER :: ZETAIT
          REAL, DIMENSION(:,:,:), POINTER :: SSZ
          REAL, DIMENSION(:), POINTER :: EPS
          REAL, DIMENSION(:), POINTER :: NUS
          REAL, DIMENSION(:), POINTER :: DELNUS
          REAL, DIMENSION(:), POINTER :: NUPLANE
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICR
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICC
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICUMCR
          REAL, DIMENSION(:,:,:,:), POINTER :: SWICUMCC
          REAL, DIMENSION(:,:,:), POINTER :: NUTOP
          REAL, DIMENSION(:,:,:), POINTER :: NUBOT
          REAL, DIMENSION(:,:,:), POINTER :: BRHS
          REAL, DIMENSION(:,:,:), POINTER :: QZEXTRA
          REAL, DIMENSION(:,:,:), POINTER :: QZEXTRACUM
          DOUBLE PRECISION, DIMENSION(:,:), POINTER :: DUM
          REAL, DIMENSION(:,:,:), POINTER :: RHSFRESH
          REAL, DIMENSION(:,:,:), POINTER :: HCOFFRESH
          INTEGER, DIMENSION(:,:,:,:), POINTER :: IPLPOS
          INTEGER, DIMENSION(:,:,:), POINTER :: IZONENR
          REAL, DIMENSION(:,:,:), POINTER ::SC1
C-----------STORAGE FOR OBSERVATION DATA
          TYPE (TSWIOBS), DIMENSION(:), POINTER :: SWIOBS
C-----------STORAGE FOR SOLVERS
          TYPE (TSWIDE4),  POINTER :: SWIDE4
          TYPE (TSWIPCG),  POINTER :: SWIPCG
          TYPE (TSWIPCGU), POINTER :: SWIPCGU
        END TYPE
        TYPE(GWFSWITYPE), SAVE:: GWFSWIDAT(10)
      END MODULE GWFSWIMODULE


      SUBROUTINE GWF2SWI1AR(In,Ibcf,Ilpf,Ihuf,Igrid)
C-----VERSION 3 26SEP2010 GWF2SWI1AR
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR SEA WATER INTRUSION PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                      LBOTM,BOTM
        USE GWFBCFMODULE,ONLY:SC1B=>SC1
        USE GWFLPFMODULE,ONLY:SC1L=>SC1
        USE GWFHUFMODULE,ONLY:SC1H=>SC1
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: In
        INTEGER, INTENT(IN) :: Ibcf
        INTEGER, INTENT(IN) :: Ilpf
        INTEGER, INTENT(IN) :: Ihuf
        INTEGER, INTENT(IN) :: Igrid
C       + + + LOCAL DEFINITIONS + + +
        CHARACTER*200 :: line
        INTEGER :: lloc, istart, istop
        INTEGER :: ierr
        INTEGER :: i, j, k, n
        INTEGER :: iz, kk
        INTEGER :: itmem
        REAL :: r
        REAL :: d
        REAL :: bbot, ttop, z
        CHARACTER*40, DIMENSION(3) :: csolver
        CHARACTER*24, DIMENSION(4) :: ANAME
        CHARACTER*24 :: ZETANAME
C       + + + DATA + + +
        DATA csolver /'                     DIRECT SOLVER (DE4)',
     2                ' PRECONDITIONED CONJUGATE GRADIENT (PCG)',
     3                'PRECONDITIONED CONJUGATE GRADIENT (PCGU)'/

        DATA ANAME(1) /'                  NUZONE'/
        DATA ANAME(2) /'                 NUPLANE'/
        DATA ANAME(3) /'                     SSZ'/
        DATA ANAME(4) /'                 IZONENR'/
C     + + + INPUT FORMATS + + +
C     + + + OUTPUT FORMATS + + +
    1   FORMAT(//1X,'SWI7 -- SWI PACKAGE, VERSION 3.1, 16-Jun-2011',
     2          /1X,45('-'),
     3          /1X,'SWI INPUT READ FROM UNIT',I3,//)
    7   FORMAT(//1X,'SWI ERROR: THE TOTAL NUMBER OF PLANES LESS THAN 1')
02200   FORMAT(//1X,'SWI DATASET 1',/,1X,52('-'),
     2    /1X,'NUMBER OF SURFACES (NPLN):                    ',1X,I5,
     3    /1X,'TOTAL NUMBER OF ZONES:                        ',1X,I5,
     4    /1X,'DENSITY DISTRIBUTION OPTION (ISTRAT):         ',1X,I5,
     5    /1X,'  VARIABLE DENSITY FLOW -- ISTRAT = 0',
     6    /1X,'  STRATIFIED FLOW       -- ISTRAT = 1',
     7    /1X,'NUMBER OF SWI OBSERVATIONS (NOBS):            ',1X,I5,
     8    /1X,'ZETA OUTPUT FILE UNIT NUMBER (ISWIZT):        ',1X,I5,
     9    /1X,'BUDGET OUTPUT FILE UNIT NUMBER (ISWICB):      ',1X,I5,
     X    /1X,'OBSERVATION OUTPUT FILE UNIT NUMBER (ISWIOBS):',1X,I5,
     1    /1X,'FRANS SCHAARS SUPER SPECIAL VARIABLE (NOPT):  ',1X,I5,
     2    /1X,52('-'))
02210   FORMAT(//1X,'SWI DATASET 2A',/,1X,52('-'),
     2    /1X,'ERROR SPECIFYING NSOLVER - VALID VALUES ARE',
     3    100(:/1X,I2,1X,'=',A40))
02220   FORMAT(//1X,'SWI DATASET 2A',/,1X,52('-'),
     2    /1X,'SOLVER (NSOLVER):                             ',1X,I5,
     3    /1X,'  SOLVER =',1X,A40,
     4    /1X,'SOLVER PRINTOUT INTERVAL (IPRSOL):            ',1X,I5,
     5    /1X,'SOLVER OUTPUT FLAG (MUTSOL):                  ',1X,I5,
     6    /1X,'  0 = PRINTING EVERY ITERATION',
     7    /1X,'  1 = LIMITED (TOTAL NUMBER OF ITERATIONS)',
     8    /1X,'  2 = NO PRINTING',
     9    /1X,'  3 = ONLY IF CONVERGENCE FAILS',
     X    /1X,52('-'))
  111   FORMAT('           ZETA PLANE ',I2)
 2100   FORMAT(//1X,'SWI PROCESS REQUIRES USE OF THE BCF, LPF,'
     2           1X,'OR HUF FLOW PACKAGES',//)
 2110   FORMAT(//1X,16X,'SWI OBSERVATION LOCATIONS',
     2          /1X,' OBSERVATION',
     3           1X,'     LAYER',1X,'       ROW',1X,'    COLUMN',
     4           1X,'      OBSNAM',
     5          /1X,58('-'))
 2120   FORMAT(1X,I12,3(1X,I10),1X,A12)
 2140   FORMAT(//1X,13X,'SWI INVALID OBSERVATION LOCATIONS',
     2          /1X,59('-'))
 2150   FORMAT(1X,' OBSERVATION',1X,I5,1X,A6,1X,I5,
     2         1X,'NOT BETWEEN',1X,I5,1X,'AND',1X,I5)
 2160   FORMAT(1X,'        NONE')
C     ------------------------------------------------------------------
C
C       + + + CODE + + +
C
C---------ALLOCATE VARIABLES - INITIALIZE IF POSSIBLE
        ALLOCATE(NPLN,ISTRAT,NZONES)
        ALLOCATE(NOBS,IOBSHEADER)
        ALLOCATE(ISWIZT,ISWICB,ISWIOBS)
        ALLOCATE(NOPT)
        ALLOCATE(NLAYSWI)
        ALLOCATE(NSOLVER,IPRSOL,MUTSOL)
        ALLOCATE(TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA)
        ALLOCATE(NADPTFLG,NADPTMX,NADPTMN,ADPTFCT)
        ALLOCATE(IADPT,IADPTMOD,ADPTVAL,SWIDELT)

        IOBSHEADER = 0
C
C---------IDENTIFY PACKAGE AND INITIALIZE
        WRITE(IOUT,1) In
C
C---------READ DATASET 1
        CALL URDCOM(In, IOUT, line)
        lloc = 1
        CALL URWORD(line, lloc, istart, istop, 2,   NPLN, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2, ISTRAT, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2,   NOBS, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2, ISWIZT, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2, ISWICB, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2,ISWIOBS, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2,   NOPT, r, IOUT, In)
C
C---------CHECK TO SEE THAT NUMBER OF PLANES IS AT LEAST 1
        IF ( NPLN.LT.1 ) THEN
          WRITE (IOUT,7)
          CALL USTOP('SWI ERROR: NUMBER OF PLANES LESS THAN 1')
        ENDIF
C
C---------CALCULATE THE TOTAL NUMBER OF ZONES
        NZONES  = NPLN + 1
        NLAYSWI = 1
C
C---------WRITE DATASET 1
        WRITE (IOUT,2200) NPLN, NZONES, ISTRAT, NOBS,
     2                    ISWIZT, ISWICB, ISWIOBS, NOPT
C
C---------READ DATASET 2A - SOLVER DATA FOR BOTH SOLVERS
        CALL URDCOM(In, IOUT, line)
        lloc = 1
        CALL URWORD(line, lloc, istart, istop, 2,NSOLVER, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2, IPRSOL, r, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2, MUTSOL, r, IOUT, In)
        IF ( NSOLVER.LT.1 .OR. NSOLVER.GT.3 ) THEN
          WRITE (IOUT,2210) (i,csolver(i),i=1,3)
          CALL USTOP('SWI ERROR: INVALID NSOLVER SPECIFIED')
        END IF
        IF ( IPRSOL.LT.1 ) IPRSOL = 999
        IF ( MUTSOL.LT.0 .OR. MUTSOL.GT.3 ) MUTSOL=0
C
C---------WRITE DATASET 2A
        WRITE (IOUT,2220) NSOLVER, csolver(NSOLVER),
     2                    IPRSOL, MUTSOL

C
C---------ALLOCATE STORAGE FOR SOLVER DATA
        ALLOCATE(SWIDE4,SWIPCG,SWIPCGU)
C
C---------READ DATASET 2B IF NECESSARY
        SELECT CASE (NSOLVER)
          CASE (2)
            CALL URDCOM(In, IOUT, line)
            lloc = 1
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCG%MXITER = i
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCG%ITER1  = i
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCG%NPCOND = i
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCG%IHCOFADD = i
            CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
            SWIPCG%ZCLOSEPCG = r
            CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
            SWIPCG%RCLOSEPCG = r
            CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
            SWIPCG%RELAXPCG = r
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCG%NBPOL = i
            CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
            SWIPCG%DAMPPCG = ABS(r)
            IF ( r.LT.0 ) THEN
              CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
              SWIPCG%DAMPPCGT = r
            ELSE
              SWIPCG%DAMPPCGT = SWIPCG%DAMPPCG
            END IF
C             WRITE DATASET 2B FOR PCG SOLVER
            WRITE (IOUT,2230) SWIPCG%MXITER, SWIPCG%ITER1,
     2                        SWIPCG%NPCOND, SWIPCG%IHCOFADD,
     3                        SWIPCG%ZCLOSEPCG,
     4                        SWIPCG%RCLOSEPCG,
     5                        SWIPCG%RELAXPCG, SWIPCG%NBPOL,
     6                        SWIPCG%DAMPPCG, SWIPCG%DAMPPCGT
          CASE (3)
            CALL URDCOM(In, IOUT, line)
            lloc = 1
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCGU%MXITER = i
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCGU%ITER1C = i
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIPCGU%NPC    = ABS(i)
            SWIPCGU%NORD   = 0
            IF ( i.EQ.-3 .OR. i.EQ.-4 ) THEN
              SWIPCGU%NORD   = 1
            END IF
            CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
            SWIPCGU%ZCLOSEPCGU = r
            CALL URWORD(line, lloc, istart, istop, 3, i, r, IOUT, In)
            SWIPCGU%RCLOSEPCGU = r
C             WRITE DATASET 2A FOR PCGU SOLVER
            WRITE (IOUT,2235) SWIPCGU%MXITER, SWIPCGU%ITER1C,
     2                        SWIPCGU%NPC, SWIPCGU%NORD,
     3                        SWIPCGU%ZCLOSEPCGU,
     4                        SWIPCGU%RCLOSEPCGU
        END SELECT
02230   FORMAT(//1X,'SWI DATASET 2B',
     2    /1X,'PRECONDITIONED CONJUGATE GRADIENT PARAMETERS',
     3    /1X,52('-'),
     4    /1X,'MAXIMUM NUMBER OF CALLS TO PCG (MXITER):      ',1X,I5,
     5    /1X,'MAXIMUM ITERATIONS PER CALL TO PCG (ITER1):   ',1X,I5,
     6    /1X,'MATRIX PRECONDITIONING TYPE (NPCOND):         ',1X,I5,
     7    /1X,'MODIFIED MATRIX ASSEMBLY OPTION (IHCOFADD):   ',1X,I5,
     8    /1X,'  1 = MODIFIED INCOMPLETE CHOLESKY',
     9    /1X,'  2 = NEUMAN POLYNOMIAL - THE MATRIX WILL BE SCALED',
     X    /1X,'ZETA CLOSURE CRITERION (ZCLOSE):    ',1X,G15.5,
     1    /1X,'RESIDUAL CLOSURE CRITERION (RCLOSE):',1X,G15.5,
     2    /1X,'RELAXATION FACTOR (RELAX):          ',1X,G15.5,
     3    /1X,'  ONLY USED WITH NPCOND = 1',
     4    /1X,'POLYNOMIAL PRECONDITIONER PARAMETER (NBPOL):  ',1X,I5,
     5    /1X,'  INTERNALLY CALCULATED IF NPCOND.NE.2',
     6    /1X,'STEADY-STATE DAMPING (DAMPPCG):     ',1X,G15.5,
     7    /1X,'TRANSIENT DAMPING (DAMPPCGT):       ',1X,G15.5,
     8    /1X,52('-'))
02235   FORMAT(//1X,'SWI DATASET 2B',
     2    /1X,'UNSTRUCTURED PRECONDITIONED ',
     3     1X,'CONJUGATE GRADIENT PARAMETERS',
     4    /1X,52('-'),
     5    /1X,'MAXIMUM NUMBER OF CALLS TO PCGU (MXITER):     ',1X,I5,
     6    /1X,'MAXIMUM ITERATIONS PER CALL TO PCGU (ITER1):  ',1X,I5,
     7    /1X,'MATRIX PRECONDITIONING TYPE (NPC):            ',1X,I5,
     8    /1X,'  1 = NO PRECONDITIONER',
     9    /1X,'  2 = JACOBI PRECONDITIONER',
     X    /1X,'  3 = ILU0 PRECONDITIONER',
     1    /1X,'  4 = MILU0 PRECONDITIONER',
     2    /1X,'  5 = NEUMAN POLYNOMIAL - THE MATRIX WILL BE SCALED',
     3    /1X,'ILU0/MILU0 PRECONDITIONER REODERING (NORD):   ',1X,I5,
     4    /1X,'  0 = NO REORDERING',
     5    /1X,'  1 = JACOBI PRECONDITIONER',
     6    /1X,'ZETA CLOSURE CRITERION (ZCLOSE):    ',1X,G15.5,
     7    /1X,'RESIDUAL CLOSURE CRITERION (RCLOSE):',1X,G15.5,
     8    /1X,52('-'))
C
C---------ALLOCATE ARRAYS
        IF ( NOBS.GT.0 ) THEN
          ALLOCATE(SWIOBS(NOBS))
        ELSE
          ALLOCATE(SWIOBS(1))
        END IF
        ALLOCATE(ZETA(NCOL,NROW,NLAY,NZONES+1))
        ALLOCATE(SSZ(NCOL,NROW,NLAY))
        ALLOCATE(IZONENR(NCOL,NROW,NLAY))
        ALLOCATE(NUS(NZONES))
        ALLOCATE(NUPLANE(NZONES+1))
        ALLOCATE(DELNUS(NZONES))
        ALLOCATE(EPS(NZONES))

        ALLOCATE(IPLPOS(NCOL,NROW,NLAY,NZONES))
        ALLOCATE(RHSFRESH(NCOL,NROW,NLAY))
        ALLOCATE(SWICR(NCOL,NROW,NLAY,NZONES))
        ALLOCATE(SWICC(NCOL,NROW,NLAY,NZONES))
        ALLOCATE(SWICUMCR(NCOL,NROW,NLAY,NZONES))
        ALLOCATE(SWICUMCC(NCOL,NROW,NLAY,NZONES))
        ALLOCATE(NUTOP(NCOL,NROW,NLAY),NUBOT(NCOL,NROW,NLAY))
        ALLOCATE(QZEXTRA(NCOL,NROW,NLAY))
        ALLOCATE(QZEXTRACUM(NCOL,NROW,NLAY))

        ALLOCATE(BRHS(NCOL,NROW,NZONES))
        ALLOCATE(HCOFFRESH(NCOL,NROW,NLAY))
        ALLOCATE(ZETAOLD(NCOL,NROW,NLAY,NZONES+1))
!        ALLOCATE(ZETAIT(NCOL,NROW,NLAY,NZONES+1))
        ALLOCATE(IBO(NCOL,NROW),SWIHCOF(NCOL,NROW))
        ALLOCATE(SWISOLVCR(NCOL,NROW))
        ALLOCATE(SWISOLVCC(NCOL,NROW))
        ALLOCATE(SWISOLVCV(NCOL,NROW))
        ALLOCATE(DUM(NCOL,NROW))
C
C-------READ DATASET 3 PARAMETERS
        CALL URDCOM(In, IOUT, line)
        lloc = 1
        CALL URWORD(line, lloc, istart, istop, 3, i, TOESLOPE, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 3, i, TIPSLOPE, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 3, i,  ZETAMIN, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 3, i,  DELZETA, IOUT, In)
        CALL URWORD(line, lloc, istart, istop, 2,  NADPTMX, r,-IOUT, In)
        IF ( NADPTMX.GT.1 ) THEN
          CALL URWORD(line, lloc, istart, istop, 2,NADPTMN, r, IOUT, In)
          IF ( NADPTMX.GT.NADPTMN ) THEN
            CALL URWORD(line,lloc,istart,istop,3,i,ADPTFCT,IOUT,In)
          ELSE 
            ADPTFCT = 1.0
          END IF
        ELSE IF ( NADPTMX.LT.1 ) THEN
          NADPTMX  = 1
          NADPTMN  = 1
          ADPTFCT = 1.0
        END IF
        IF ( NADPTMN.LT.1 ) THEN
          NADPTMN = 1
        END IF
        NADPTFLG = 0
        IF ( NADPTMX.GT.NADPTMN ) NADPTFLG = 1

        IADPT      = NADPTMN
        IADPTMOD   = -1
        ADPTVAL    = 1.0

        WRITE (IOUT,2240) TOESLOPE,TIPSLOPE,ZETAMIN,DELZETA,
     2                    NADPTMN,NADPTMX
        IF ( NADPTFLG.NE.0 ) THEN
          WRITE (IOUT,2245) ADPTFCT
        END IF
        WRITE (IOUT,2250)
        
02240   FORMAT(//1X,'SWI DATASET 3',/,1X,52('-'),
     2    /1X,'MAXIMUM TOE SURFACE SLOPE (TOESLOPE):    ',1X,G10.3,
     3    /1X,'MAXIMUM TIP SURFACE SLOPE (TIPSLOPE):    ',1X,G10.3,
     4    /1X,'MINIMUM ZETA SURFACE (ZETAMIN):          ',1X,G10.3,
     5    /1X,'ZETA INCREMENT (DELZETA):                ',1X,G10.3,
     6    /1X,'MINIMUM NUMBER OF SUB-TIME STEPS (NADPTMN):   ',1X,I5,
     7    /1X,'MAXIMUM NUMBER OF SUB-TIME STEPS (NADPTMX):   ',1X,I5)
02245   FORMAT(
     2    /1X,'MAXIMUM ZETA SLOPE FRACTION (ADPTFCT):   ',1X,G10.4)
02250   FORMAT(1X,52('-'))

C---------READ DENSITY DATA BASED ON USER-SPECIFIED DENSITY MODEL
C         CONSTANT ZONE DENSITY MODEL
        DENSITYMODEL: IF (ISTRAT.EQ.1) THEN
C
C-----------READ NU FOR EACH ZONE (NZONES)
          CALL U1DREL(NUS,ANAME(1),NZONES,In,IOUT)
C-----------SET EPS EQUAL TO ZERO
          DO iz = 1, NZONES
            EPS(iz) = 0.0
          END DO
C         LINEAR ZONE DENSITY MODEL
        ELSEIF (ISTRAT.EQ.0) THEN
C
C-----------READ NUPLANE FOR EACH PLANE (NZONES+1)
          CALL U1DREL(NUPLANE,ANAME(2),NZONES+1,In,IOUT)
C
C-----------CALCULATE NU AND EPS FROM NUPLANE
          DO iz = 1, NZONES
            NUS(iz) = 0.5*(NUPLANE(iz)+NUPLANE(iz+1))
            EPS(iz) = (NUPLANE(iz+1)-NUPLANE(iz)) / 6
          END DO
        END IF DENSITYMODEL
C
C---------CALCULATE DELNUS FROM NUS
        DELNUS(1)=NUS(1);
        DO iz = 2, NZONES
          DELNUS(iz)=(NUS(iz)-NUS(iz-1))
        END DO
C
C---------READ ZETA FOR EACH PLANE
        IZ_ZETA: DO iz=2,NZONES
          K_ZETA: DO k=1,NLAY
            kk = k
            WRITE(ZETANAME,111) iz-1
            CALL U2DREL(ZETA(1,1,k,iz),ZETANAME,NROW,NCOL,kk,In,IOUT)
C
C-------------LOOP THROUGH EACH ROW AND COLUMN
C             RESET ZETA IF SPECIFIED ZETA IS GREATER THAN THE TOP
C             OF THE CURRENT CELL OR LESS THAN THE BOTTOM OF THE CURRENT
C             CELL - IF ZETA IS RESET TO THE TOP OR BOTTOM BASED ON IF
C             THE SPECIFIED ZETA VALUE IS CLOSER TO THE TOP OR BOTTOM
            d = 0.001
            I_ZETA: DO i = 1, NROW
              J_ZETA: DO j = 1, NCOL
                bbot = BOTM(j,i,LBOTM(k))
                ttop = BOTM(j,i,LBOTM(k)-1)
                z = ZETA(j,i,k,iz)
                IF (z.NE.bbot) THEN
                  IF (z.LT.(bbot+d)) THEN
                    ZETA(j,i,k,iz) = bbot
                  END IF
                END IF
                IF (z.GT.(ttop-d)) ZETA(j,i,k,iz) = ttop
                ZETAOLD(j,i,k,iz) = ZETA(j,i,k,iz)
!                ZETAIT(j,i,k,iz)  = ZETA(j,i,k,iz)
              END DO J_ZETA
            END DO I_ZETA

          END DO K_ZETA
        END DO IZ_ZETA
C
C---------READ SSZ FOR EACH LAYER
        K_SSZ: DO k = 1, NLAY
          kk=k
          CALL U2DREL(SSZ(1,1,k),ANAME(3),NROW,NCOL,kk,In,IOUT)
        END DO K_SSZ
C
C---------READ IZONENR FOR EACH LAYER
        K_IZONENR: DO k = 1, NLAY
          kk=k
          CALL U2DINT(IZONENR(1,1,kk),ANAME(4),NROW,NCOL,kk,In,IOUT)
        END DO K_IZONENR
C
C---------READ OBSERVATION DATA
        IF ( NOBS.GT.0 ) THEN
          DO n = 1, NOBS
C-------------READ OBSERVATION DATA
            CALL URDCOM(In, IOUT, line)
            lloc = 1
            CALL URWORD(line, lloc, istart, istop, 0, i, r, IOUT, In)
            SWIOBS(n)%OBSNAM = line(istart:istop)
            CALL URWORD(line, lloc, istart, istop, 2, k, r, IOUT, In)
            SWIOBS(n)%KLAY = k
            CALL URWORD(line, lloc, istart, istop, 2, i, r, IOUT, In)
            SWIOBS(n)%IROW = i
            CALL URWORD(line, lloc, istart, istop, 2, j, r, IOUT, In)
            SWIOBS(n)%JCOL = j
          END DO
C-------------WRITE OBSERVATION DATA TO IOUT
          WRITE (IOUT,2110)
          DO n = 1, NOBS
            WRITE (IOUT,2120)
     2        n, SWIOBS(n)%KLAY, SWIOBS(n)%IROW, SWIOBS(n)%JCOL,
     3        ADJUSTR(TRIM(SWIOBS(n)%OBSNAM))
          END DO
C-----------CHECK FOR INVALID OBSERVATION LOCATIONS
          ierr = 0
          WRITE(IOUT,2140)
          DO n = 1, NOBS
            k = SWIOBS(n)%KLAY
            i = SWIOBS(n)%IROW
            j = SWIOBS(n)%JCOL
            IF ( k.LT.1 .OR. k.GT.NLAY ) THEN
              ierr = ierr + 1
              WRITE (IOUT,2150) n, 'LAYER ', k, 1, NLAY
            END IF
            IF ( i.LT.1 .OR. i.GT.NROW ) THEN
              ierr = ierr + 1
              WRITE (IOUT,2150) n, 'ROW   ', i, 1, NROW
            END IF
            IF ( j.LT.1 .OR. j.GT.NCOL ) THEN
              ierr = ierr + 1
              WRITE (IOUT,2150) n, 'COLUMN', j, 1, NCOL
            END IF
          END DO
          IF ( ierr.GT.0 ) THEN
            CALL USTOP('SWI ERROR: INVALID OBSERVATION LOCATIONS')
          ELSE
            WRITE (IOUT,2160)
          END IF
        END IF
C
C---------POINT SWI SC1 TO APPROPRIATE STORAGE PARAMETER
        IF ( Ibcf.GT.0 ) THEN
          SC1 => SC1B
        ELSE IF ( Ilpf.GT.0 ) THEN
          SC1 => SC1L
        ELSE IF ( Ihuf.GT.0 ) THEN
          SC1 => SC1H
        ELSE
          WRITE (IOUT,2100)
          CALL USTOP('FLOW PACKAGE SPECIFIED INCONSISTENT WITH SWI')
        END IF
C
C---------ALLOCATE SPACE FOR THE SOLVER ARRAYS
C         DIRECT SOLVER (DE4)
        SOLVERPARAM: IF ( NSOLVER.EQ.1 ) THEN
          SWIDE4%ITMX = 1
!-----------SPECIFY DEFAULT PARAMETERS
          SWIDE4%MXITER = SWIDE4%ITMX
          SWIDE4%NITERDE4 = 1
          SWIDE4%ACCLDE4 = 1.0
          SWIDE4%HCLOSEDE4 = 1.0E-05
          SWIDE4%IFREQ = 3
          SWIDE4%DELTL = 1.
C           INITIALIZE DIMENSION PARAMETERS
          SWIDE4%NODES=NCOL*NROW*NLAYSWI
          SWIDE4%NHALFU=(SWIDE4%NODES-1)/2 + 1
          SWIDE4%NHALFL=SWIDE4%NODES-SWIDE4%NHALFU
          SWIDE4%NBWL = 0
          SWIDE4%NUPL = 0
          SWIDE4%NLOWL = 0
          SWIDE4%ID4DIM = 7
C           CALCULATE SOLVER DIMENSIONS
          IF(NLAYSWI.LE.NCOL .AND. NLAYSWI.LE.NROW) THEN
             IF(NLAYSWI.EQ.1) SWIDE4%ID4DIM=5
             IF(NCOL.GE.NROW) THEN
                SWIDE4%ID4DIR=1
                SWIDE4%NBWGRD=NROW*NLAYSWI+1
             ELSE
                SWIDE4%ID4DIR=2
                SWIDE4%NBWGRD=NCOL*NLAYSWI+1
             END IF
          ELSE IF(NROW.LE.NCOL .AND. NROW.LE.NLAYSWI) THEN
             IF(NROW.EQ.1) SWIDE4%ID4DIM=5
             IF(NCOL.GE.NLAYSWI) THEN
                SWIDE4%ID4DIR=3
                SWIDE4%NBWGRD=NROW*NLAYSWI+1
             ELSE
                SWIDE4%ID4DIR=4
                SWIDE4%NBWGRD=NROW*NCOL+1
             END IF
          ELSE
             IF(NCOL.EQ.1) SWIDE4%ID4DIM=5
             IF(NROW.GE.NLAYSWI) THEN
                SWIDE4%ID4DIR=5
                SWIDE4%NBWGRD=NCOL*NLAYSWI+1
             ELSE
                SWIDE4%ID4DIR=6
                SWIDE4%NBWGRD=NCOL*NROW+1
             END IF
          END IF
          SWIDE4%MXUP=SWIDE4%NHALFU
          SWIDE4%MXLOW=SWIDE4%NHALFL
          SWIDE4%MXBW=SWIDE4%NBWGRD + 4
          SWIDE4%MXEQ=SWIDE4%MXUP+SWIDE4%MXLOW
C         ALLOCATE DE4 ARRAYS
          ALLOCATE (SWIDE4%AU(SWIDE4%ID4DIM,SWIDE4%MXUP))
          ALLOCATE (SWIDE4%IUPPNT(SWIDE4%ID4DIM,SWIDE4%MXUP))
          ALLOCATE (SWIDE4%AL(SWIDE4%MXBW,SWIDE4%MXLOW))
          ALLOCATE (SWIDE4%IEQPNT(NCOL,NROW,NLAYSWI))
          ALLOCATE (SWIDE4%D4B(SWIDE4%MXEQ))
          ALLOCATE (SWIDE4%LRCHDE4(3,SWIDE4%ITMX))
          ALLOCATE (SWIDE4%HDCGDE4(SWIDE4%ITMX))
C           INITIALIZE FLOAT ARRAYS
          DO i = 1, SWIDE4%MXUP
            DO j = 1, SWIDE4%ID4DIM
              SWIDE4%AU(j,i) = 0.0
            END DO
          END DO
          DO i = 1, SWIDE4%MXLOW
            DO j = 1, SWIDE4%MXBW
              SWIDE4%AL(j,i) = 0.0
            END DO
          END DO
          DO i = 1, SWIDE4%MXEQ 
            SWIDE4%D4B(i) = 0.0
          END DO
C         PRECONDITIONED CONJUGATE GRADIENT SOLVER (PCG4)
        ELSE IF ( NSOLVER.EQ.2 ) THEN
C           INITIALIZE DIMENSION PARAMETERS
          SWIPCG%NODES=NCOL*NROW*NLAYSWI
C---------ALLOCATE SPACE FOR THE PCG ARRAYS
          ALLOCATE (SWIPCG%VPCG(NCOL,NROW,NLAYSWI))
          ALLOCATE (SWIPCG%SS(NCOL,NROW,NLAYSWI))
          ALLOCATE (SWIPCG%P(NCOL,NROW,NLAYSWI))
          ALLOCATE (SWIPCG%HPCG(NCOL,NROW,NLAYSWI))
          ALLOCATE (SWIPCG%CD(NCOL,NROW,NLAYSWI))
          IF(SWIPCG%NPCOND.EQ.2) THEN
             ALLOCATE (SWIPCG%HCSV(NCOL,NROW,NLAYSWI))
          ELSE
             ALLOCATE (SWIPCG%HCSV(1,1,1))
          END IF
          itmem=SWIPCG%MXITER*SWIPCG%ITER1
          ALLOCATE (SWIPCG%HCHG(itmem))
          ALLOCATE (SWIPCG%LHCH(3,itmem))
          ALLOCATE (SWIPCG%RCHG(itmem))
          ALLOCATE (SWIPCG%LRCHPCG(3,itmem))
          ALLOCATE (SWIPCG%IT1(itmem))
C           INITIALIZE FLOAT ARRAYS - FIX -- WILL CAUSE STACK ERRORS
          DO k = 1, NLAYSWI
            DO i = 1, NROW
              DO j = 1, NCOL
                SWIPCG%VPCG(j,i,k) = 0.0
                SWIPCG%SS(j,i,k)   = 0.0
                SWIPCG%P(j,i,k)    = 0.0
                SWIPCG%HPCG(j,i,k) = 0.0
                SWIPCG%CD(j,i,k)   = 0.0
                IF(SWIPCG%NPCOND.EQ.2) THEN
                  SWIPCG%HCSV(j,i,k) = 0.0
                END IF
              END DO
            END DO
          END DO
C         UNSTRUCTURED PRECONDITIONED CONJUGATE GRADIENT SOLVER (PCGU)
        ELSE  IF ( NSOLVER.EQ.3 ) THEN
C
C-----------ALLOCATE SPACE FOR SOLVER CONVERGENCE ARRAYS
          itmem = SWIPCGU%MXITER * SWIPCGU%ITER1C
          ALLOCATE (SWIPCGU%HCHGPCGU(itmem))
          ALLOCATE (SWIPCGU%LHCHPCGU(3,itmem))
          ALLOCATE (SWIPCGU%RCHGPCGU(itmem))
          ALLOCATE (SWIPCGU%LRCHPCGU(3,itmem))
          ALLOCATE (SWIPCGU%IT1PCGU(itmem))
C
C-----------CALCULATE NUMBER OF ENTRIES IN MODEL GRID
          SWIPCGU%NNZC  = 0
          SWIPCGU%NIAC  = 0
          SWIPCGU%NIAPC = 0
          SWIPCGU%NIWC  = 0
          SWIPCGU%NODES=NCOL*NROW*NLAYSWI
          ALLOCATE( SWIPCGU%NODEC(NCOL,NROW,NLAYSWI) )
          CALL PCGU7CCN(NLAYSWI,NROW,NCOL,
     2                  SWIPCGU%NIAC,SWIPCGU%NNZC,SWIPCGU%NODEC)
C
C-----------ALLOCATE AND INITIALIZE COMPRESSED ROW STORAGE VECTORS
          ALLOCATE(SWIPCGU%AC(SWIPCGU%NNZC))
          SWIPCGU%NIAPC = SWIPCGU%NNZC
          IF ( SWIPCGU%NPC.EQ.1 ) THEN
            SWIPCGU%NIAPC = 1
          ELSE IF ( SWIPCGU%NPC.EQ.2 ) THEN
            SWIPCGU%NIAPC = SWIPCGU%NIAC
          END IF
          ALLOCATE(SWIPCGU%TAPC(SWIPCGU%NIAPC))
          ALLOCATE(SWIPCGU%APC(SWIPCGU%NIAPC))
          ALLOCATE(SWIPCGU%IAC(SWIPCGU%NIAC+1))
          ALLOCATE(SWIPCGU%JAC(SWIPCGU%NNZC))
          ALLOCATE(SWIPCGU%IUC(SWIPCGU%NIAC))
          ALLOCATE(SWIPCGU%IXMAP(SWIPCGU%NIAC))
C           ALLOCATE SPACE FOR ILU0 AND MILU0 NON-ZERO ROW ENTRY VECTOR
          IF ( SWIPCGU%NPC.EQ.3 .OR.
     2         SWIPCGU%NPC.EQ.4 .OR.
     3         SWIPCGU%NPC.EQ.5 ) THEN
            SWIPCGU%NIWC = SWIPCGU%NIAC
          ELSE
            SWIPCGU%NIWC = 1
          END IF
          ALLOCATE(SWIPCGU%IWC(SWIPCGU%NIWC))
C           ALLOCATE SPACE FOR REODERING ILU0 AND MILU0 NON-ZERO ELEMENTS
          SWIPCGU%NIAPROC  = SWIPCGU%NIAC
          SWIPCGU%NNZAPROC = SWIPCGU%NNZC
          IF ( SWIPCGU%NORD.LT.1 ) THEN
            SWIPCGU%NIAPROC  = 1
            SWIPCGU%NNZAPROC = 1
          END IF
          ALLOCATE(SWIPCGU%IAPC(SWIPCGU%NIAPROC+1))
          ALLOCATE(SWIPCGU%JAPC(SWIPCGU%NNZAPROC))
          ALLOCATE(SWIPCGU%IRAPC(SWIPCGU%NNZAPROC))
          ALLOCATE(SWIPCGU%IAPCMAP(SWIPCGU%NNZAPROC))
C           ALLOCATE WORKING VECTORS FOR PCGU SOLVER
          ALLOCATE(SWIPCGU%BC(SWIPCGU%NIAC),SWIPCGU%XC(SWIPCGU%NIAC))
          ALLOCATE(SWIPCGU%DC(SWIPCGU%NIAC),SWIPCGU%PC(SWIPCGU%NIAC))
          ALLOCATE(SWIPCGU%QC(SWIPCGU%NIAC),SWIPCGU%ZC(SWIPCGU%NIAC))
          ALLOCATE(SWIPCGU%SCL(SWIPCGU%NIAC))
C           INITIALIZE PCGU SOLVER VECTORS
          DO n = 1, SWIPCGU%NNZC
            SWIPCGU%AC(n)  = 0.0D0
            SWIPCGU%JAC(n) = 0
          END DO
          DO n = 1, SWIPCGU%NIAPC
            SWIPCGU%APC(n) = 0.0D0
          END DO
          DO n = 1, SWIPCGU%NIAC+1
            SWIPCGU%IAC(n) = 0
          END DO
          DO n = 1, SWIPCGU%NIAC
            SWIPCGU%IUC(n)   = 0
            SWIPCGU%IXMAP(n) = 0
            SWIPCGU%BC(n)    = 0.0D0
            SWIPCGU%XC(n)    = 0.0D0
C             WORKING ARRAYS
            SWIPCGU%DC(n)    = 0.0D0
            SWIPCGU%PC(n)    = 0.0D0
            SWIPCGU%QC(n)    = 0.0D0
            SWIPCGU%ZC(n)    = 0.0D0
            SWIPCGU%SCL(n)   = 0.0D0
          END DO
          DO n = 1, SWIPCGU%NIWC
            SWIPCGU%IWC(n)   = 0
          END DO
C
C-----------FILL IA AND JA
          CALL PCGU7FCN(SWIPCGU%NORD,NLAYSWI,NROW,NCOL,
     2      SWIPCGU%NIAC,SWIPCGU%NNZC,SWIPCGU%NLAPC,SWIPCGU%NUAPC,
     2      SWIPCGU%NIAPROC,SWIPCGU%NNZAPROC,SWIPCGU%NODEC,
     3      SWIPCGU%IAC,SWIPCGU%JAC,SWIPCGU%IUC,SWIPCGU%IXMAP,
     4      SWIPCGU%IAPC,SWIPCGU%JAPC,SWIPCGU%IRAPC,SWIPCGU%IAPCMAP)
        END IF SOLVERPARAM
C
C---------SET POINTERS FOR GRID
        CALL SGWF2SWI1PSV(Igrid)
C
C---------RETURN
        RETURN
      END SUBROUTINE GWF2SWI1AR

C
C-------
      SUBROUTINE GWFSWI1HC(Igrid)
C-----VERSION 3 26SEP2010 GWFSWI1HC
C     ******************************************************************
C     COPY HCOF TO SWI HCOFFRESH PRIOR TO ADDING MODFLOW BOUNDARY
C     CONDITION TERMS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:NCOL,NROW,NLAY,HCOF
        USE GWFSWIMODULE,ONLY:HCOFFRESH
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Igrid
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k
C     ------------------------------------------------------------------
C
C---------SET SWI POINTERS FOR GRID
        CALL SGWF2SWI1PNT(Igrid)
C
C---------HCOFFRESH TO HCOF - USED TO CALCULATE BOUNDARY FLUXES
        DO k = 1, NLAY
          DO i = 1, NROW
            DO j = 1, NCOL
              HCOFFRESH(j,i,k) = HCOF(j,i,k)
            END DO
          END DO
        END DO
C
C---------RETURN
        RETURN
      END SUBROUTINE GWFSWI1HC

      SUBROUTINE GWF2SWI1AD(Kper,Igrid)
C     ******************************************************************
C     COMPUTE HEAD FOR TIME STEP AT EACH TIME-VARIANT SPECIFIED HEAD
C     CELL.
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Kper
        INTEGER, INTENT(IN) :: Igrid
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k, n
        INTEGER :: iz
C     ------------------------------------------------------------------
C       + + + CODE + + +
C
C---------SET POINTERS FOR GRID
        CALL SGWF2SWI1PNT(Igrid)
C
C----------COPY ZETA TO ZETAOLD
        DO k = 1, NLAY
          DO i = 1, NROW
            DO j = 1, NCOL
              DO iz = 2, NZONES
                ZETAOLD(j,i,k,iz) = ZETA(j,i,k,iz)
!                ZETAIT(j,i,k,iz)  = ZETA(j,i,k,iz)
              END DO
            END DO
          END DO
        END DO
C
C---------RETURN
        RETURN
      END SUBROUTINE GWF2SWI1AD

      SUBROUTINE GWF2SWI1FM(Kkstp,Kkper,Kkiter,Igrid)
C-----VERSION 3 26SEP2010 GWF2SWI1FM
C     ******************************************************************
C     ADD SWI TERMS TO RHS AND HCOF
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                        LBOTM,BOTM,
     3                        CR,CC,CV,HCOF,RHS,
     4                        DELR,DELC,IBOUND,HNEW
        USE GWFBASMODULE, ONLY: DELT,HDRY,TOTIM
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Kkstp
        INTEGER, INTENT(IN) :: Kkper
        INTEGER, INTENT(IN) :: Kkiter
        INTEGER, INTENT(IN) :: Igrid
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: ierr
        INTEGER :: i, j, k
        INTEGER :: iz
        INTEGER :: iswi
        INTEGER :: izrev
        REAL :: bbot, ttop, hhd
        REAL :: sumthickrf, sumthickff
        REAL :: dzeta1, dzeta2
        REAL :: thickrf, thickff
        REAL :: sumswicr, sumswicc
        REAL :: headdiff
        REAL :: switfact
        REAL :: rhoc
        REAL :: r, f
C     ------------------------------------------------------------------
C       + + + CODE + + +
C
C---------SET POINTERS FOR GRID
        CALL SGWF2SWI1PNT(Igrid)
C
C--------SET FIRST AND LAST PLANE EQUAL TO TOP AND BOTTOM OF THE AQUIFER
        ierr = 0
        DO k = 1, NLAY
          DO i = 1, NROW
            DO j = 1, NCOL
              bbot=BOTM(j,i,LBOTM(k))
              ttop=BOTM(j,i,LBOTM(k)-1)
              IF (bbot.GT.ttop) THEN
                WRITE(IOUT,35) k,i,j
                WRITE(IOUT,36) ttop,bbot
                ierr = ierr + 1
              END IF
!              IF ((LAYCON(k).NE.0).OR.(LAYCON(k).NE.0)) THEN
!                HHD=HNEW(j,i,k)
!                IF(HHD.LT.ttop) ttop=HHD
!              END IF
C              hhd=HNEW(j,i,k)
C              IF(hhd.LT.ttop) ttop=hhd
              ZETA(j,i,k,1)=ttop
              ZETA(j,i,k,NZONES+1)=bbot
!              ZETAIT(j,i,k,1)=ttop
!              ZETAIT(j,i,k,NZONES+1)=bbot
              ZETAOLD(j,i,k,1)=ttop
              ZETAOLD(j,i,k,NZONES+1)=bbot
            END DO
          END DO
        END DO
   35   FORMAT(1X,/1X,'Negative cell thickness at (layer,row,col)',
     1         I4,',',I4,',',I4)
   36   FORMAT(1X,'Top elevation, bottom elevation:',1P,2G13.5)
C---------TERMINATE IF ANY CELL THICKNESSES ARE NEGATIVE
        IF (ierr.GT.0 ) THEN
          CALL USTOP('ERROR SWI: NEGATIVE CELL THICKNESSES')
        END IF
C
C-------SET ZETA TO ZETA VALUE AT THE END OF THE LAST SWI TIME STEP IN THE
C       LAST MODFLOW TIME STEP (ZETAOLD). ZETAOLD IS ASSGNED IN GWF2SWI1AD
        DO k = 1, NLAY
          DO i = 1, NROW
            DO j = 1, NCOL
              QZEXTRACUM(j,i,k) = 0.0
              DO iz = 2, NZONES
!                ZETA(j,i,k,iz) = ZETAIT(j,i,k,iz)
                ZETA(j,i,k,iz) = ZETAOLD(j,i,k,iz)
              END DO
            END DO
          END DO
        END DO
C
C---------COPY RHS TO RHSFRESH FOR CALCULATING BOUNDARY FLUXES WHEN UPDATING ZETA
        KRHS: DO k=1,NLAY
          IRHS: DO i=1,NROW
            JRHS: DO j=1,NCOL
              RHSFRESH(j,i,k) = RHS(j,i,k)
            END DO JRHS
          END DO IRHS
        END DO KRHS
C
C-------EVALUATE IF THE NUMBER OF SWI TIME STEPS PER MODFLOW
C       TIME STEP SHOULD BE REDUCED. REDUCTION IN THE NUMBER
C       OF SWI TIME STEPS PER MODFLOW TIME STEP CAN ONLY 
C       OCCUR ON THE FIRST CALL TO GWF2SWI1FM IN EACH MODFLOW
C       TIME STEP
        switfact = 1.0 / REAL( IADPT, 4 )
        SWIDELT  = DELT * switfact
!        switfact = 1.0
!        SWIDELT  = DELT
        IF ( NADPTFLG.NE.0 ) THEN
          IF ( IADPT.GT.1 .AND. Kkiter.EQ.1 ) THEN
            r = REAL( IADPT, 4 )
            IF ( r.GT.NADPTMN .AND. IADPTMOD.LT.0 ) THEN
              r = r / 2.0
              IF ( r.LT.NADPTMN ) THEN
                r = NADPTMN
              END IF
              IADPT    = CEILING( r )
              switfact = 1.0 / REAL( IADPT, 4 )
              SWIDELT  = DELT * switfact
              WRITE (IOUT,2000) IADPT, DELT, SWIDELT
            END IF
C-------------RESET THE FLAG FOR REDUCING THE NUMBER OF
C             SWI TIME STEPS (IADPTMOD)
            IADPTMOD = -1
          END IF
        END IF
2000    FORMAT(1X,'IADPT      :',1X,I10,
     2         1X,'DELT       :',1X,G10.3,
     3         1X,'SWI DELT   :',1X,G10.3)
C
C-------ADAPTIVE TIME STEP LOOP
        iswi = 1
        ADAPTIVE: DO
C
C---------COMPUTE HORIZONTAL (CUMULATIVE) CONDUCTANCE FOR EACH ZONE
        KCALC_CCCR: DO k = 1, NLAY
          ICALC_CCCR: DO i = 1, NROW
            JCALC_CCCR: DO j = 1, NCOL
              sumthickrf = 0.
              sumthickff = 0.
              IZCALC_CCCR: DO iz = 1, NZONES
                IF(j.NE.NCOL) THEN
                  dzeta1 = ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)
                  dzeta2 = ZETA(j+1,i,k,iz)-ZETA(j+1,i,k,iz+1)
                  IF (dzeta1.LE.0 .OR. dzeta2.LE.0) THEN
                    thickrf = 0.
                  ELSE
                    thickrf = (dzeta1*DELR(j+1)+dzeta2*DELR(j))/
     &                        (DELR(j)+DELR(j+1))
                  END IF
                  SWICR(j,i,k,iz) = CR(j,i,k)*thickrf
                  sumthickrf = sumthickrf+thickrf
                ELSE
                  SWICR(j,i,k,iz) = 0.
                  sumthickrf = 1.
                END IF
                IF(i.NE.NROW) THEN
                  dzeta1 = ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)
                  dzeta2 = ZETA(j,i+1,k,iz)-ZETA(j,i+1,k,iz+1)
                    IF (dzeta1.LE.0 .OR. dzeta2.LE.0) THEN
                      thickff = 0.
                    ELSE
                      thickff = (dzeta1*DELC(i+1)+dzeta2*DELC(i))/
     &                          (DELC(i)+DELC(i+1))
                    END IF
                    SWICC(j,i,k,iz) = CC(j,i,k)*thickff
                    sumthickff = sumthickff + thickff
                ELSE
                  SWICC(j,i,k,iz) = 0.
                  sumthickff = 1.
                END IF
              END DO IZCALC_CCCR
C
C...............FLEVO
              DO iz=1,NZONES
                IF (sumthickrf.NE.0) THEN
C                  sumthickrf=1
                  SWICR(j,i,k,iz) = SWICR(j,i,k,iz)/sumthickrf
                END IF
                IF (sumthickff.NE.0) THEN
C                  sumthickrf=1
                  SWICC(j,i,k,iz) = SWICC(j,i,k,iz)/sumthickff
                END IF
              END DO
C...............FLEVO
C
              sumswicr = 0.
              sumswicc = 0.
              DO iz=1,NZONES
                sumswicr = sumswicr + SWICR(j,i,k,NZONES+1-iz)
                sumswicc = sumswicc + SWICC(j,i,k,NZONES+1-iz)
                SWICUMCR(j,i,k,NZONES+1-iz) = sumswicr
                SWICUMCC(j,i,k,NZONES+1-iz) = sumswicc
              END DO
            END DO JCALC_CCCR
          END DO ICALC_CCCR
        END DO KCALC_CCCR
C
C---------SET FLAG FOR LOCATION OF THE ZETA SURFACE RELATIVE
C         TO THE TOP AND BOTTOM OF A CELL
C         IPLPOS=1 AT TOP, IPLPOS=2 AT BOTTOM, IPLPOS=0 IN BETWEEN
        K_IPLPOS: DO k=1,NLAY
          I_IPLPOS: DO i=1,NROW
            J_IPLPOS: DO j=1,NCOL
              IPLPOS(j,i,k,1) = 0
              IZ_IPLPOS: DO iz=2,NZONES
                IF (ZETA(j,i,k,iz).GE.ZETA(j,i,k,1)) THEN
                  IPLPOS(j,i,k,iz) = 1
                ELSEIF (ZETA(j,i,k,iz).LE.ZETA(j,i,k,NZONES+1)) THEN
                  IPLPOS(j,i,k,iz) = 2
                ELSE
                  IPLPOS(j,i,k,iz) = 0
                END IF
              END DO IZ_IPLPOS
            END DO J_IPLPOS
          END DO I_IPLPOS
        END DO K_IPLPOS
C
C---------SET ZONE CONDUCTANCES TO ZERO IF NOT CELL ON EITHER SIDE HAS ACTIVE ZONE
        KACTIVE: DO k=1,NLAY
          IACTIVE: DO i=1,NROW
            JACTIVE: DO j=1,NCOL
              IZACTIVE: DO iz = 2, NZONES-1
                IF (j.NE.NCOL) THEN
                  IF ((IPLPOS(j,i,k,iz)+IPLPOS(j+1,i,k,iz)+
     2                 IPLPOS(j,i,k,iz+1)+
     3                 IPLPOS(j+1,i,k,iz+1)).NE.0) THEN
                    SWICR(j,i,k,iz) = 0.
                  END IF
                END IF
                IF (i.NE.NROW) THEN
                  IF ((IPLPOS(j,i,k,iz)+IPLPOS(j,i+1,k,iz)+
     2                 IPLPOS(j,i,k,iz+1)+
     3                 IPLPOS(j,i+1,k,iz+1)).NE.0) THEN
                    SWICC(j,i,k,iz) = 0.
                  END IF
                END IF
              END DO IZACTIVE
            END DO JACTIVE
          END DO IACTIVE
        END DO KACTIVE
C
C---------COMPUTE NUTOP AND NUBOT
        KCALCNU: DO k=1,NLAY
          ICALCNU: DO i=1,NROW
            JCALCNU: DO j=1,NCOL
              NUTOP(j,i,k) = NUS(1)
              NUBOT(j,i,k) = NUS(NZONES)
              IZCALCNU: DO iz=2,NZONES
                IF (IPLPOS(j,i,k,iz).EQ.1) THEN
                  NUTOP(j,i,k) = NUTOP(j,i,k)+DELNUS(iz)
                END IF
                izrev = NZONES - iz+2
                IF (IPLPOS(j,i,k,izrev).EQ.2) THEN
                  NUBOT(j,i,k)=NUBOT(j,i,k)-DELNUS(izrev)
                END IF
              END DO IZCALCNU
            END DO JCALCNU
          END DO ICALCNU
        END DO KCALCNU
C
C---------COMPUTE QZEXTRA
        KQZE: DO k = 1, NLAY
          IQZE: DO i = 1, NROW
            JQZE: DO j = 1, NCOL
              QZEXTRA(j,i,k) = 0.0
              headdiff = 0.
              IF (k.NE.1) THEN
                IZQZE: DO iz = 2, NZONES+1
                  headdiff = headdiff -
     &              NUS(iz-1) * (ZETA(j,i,k-1,iz)-ZETA(j,i,k-1,iz-1))
                END DO IZQZE
                QZEXTRA(j,i,k) = -CV(j,i,k-1) *
     &            (headdiff+0.5*(ZETA(j,i,k-1,NZONES+1)-ZETA(j,i,k,1)) *
     &            (NUBOT(j,i,k-1)+NUTOP(j,i,k)))
              END IF
            END DO JQZE
          END DO IQZE
        END DO KQZE
C
C-------ADD QZEXTRA TO RHS
        KQZEXTRA: DO k=1,NLAY
          IQZEXTRA: DO i=1,NROW
            JQZEXTRA: DO j=1,NCOL
              !EERST=RHS(j,i,k)
              RHS(j,i,k) = RHS(j,i,k) + 
     &            QZEXTRA(j,i,k) * switfact 
              IF (k.NE.NLAY) THEN
                RHS(j,i,k) = RHS(j,i,k) - 
     &            QZEXTRA(j,i,k+1) * switfact 
              END IF
C---------------ADD QZEXTRA FOR CURRENT SWI TIME STEP TO QZEXTRACUM
              QZEXTRACUM(j,i,k) = QZEXTRACUM(j,i,k) + 
     2                            switfact * QZEXTRA(j,i,k)
              !LATER=RHS(j,i,k)
            END DO JQZEXTRA
          END DO IQZEXTRA
        END DO KQZEXTRA

C
C---------ADD TO RHS DUE TO DENSITY DIFFERENCES
        KDDIFF: DO k=1,NLAY
          IDDIFF: DO i=1,NROW
            JDDIFF: DO j=1,NCOL
              !EERST=RHS(j,i,k)
              IZDDIFF: DO iz=1,NZONES
C---------------LEFT FACE
              IF (j.NE.1) THEN
                rhoc = DELNUS(iz) * SWICUMCR(j-1,i,k,iz) *
     &            ( ZETA(j-1,i,k,iz) - ZETA(j,i,k,iz) ) +
     &            EPS(iz) * SWICR(j-1,i,k,iz) *
     &            ( ( ZETA(j-1,i,k,iz) - ZETA(j-1,i,k,iz+1) ) -
     &              ( ZETA(j,i,k,iz)   - ZETA(j,i,k,iz+1)   ) )
                RHS(j,i,k) = RHS(j,i,k) - rhoc * switfact 
              END IF
C---------------RIGHT FACE
              IF (j.NE.NCOL) THEN
                rhoc = DELNUS(iz) * SWICUMCR(j,i,k,iz) *
     &            ( ZETA(j+1,i,k,iz) - ZETA(j,i,k,iz) ) +
     &            EPS(iz) * SWICR(j,i,k,iz) *
     &            ( ( ZETA(j+1,i,k,iz) - ZETA(j+1,i,k,iz+1) ) -
     &              ( ZETA(j,i,k,iz)   - ZETA(j,i,k,iz+1)   ) )
                RHS(j,i,k) = RHS(j,i,k) - rhoc * switfact  
              END IF
C---------------BACK FACE--------------------
              IF (i.NE.1) THEN
                rhoc = DELNUS(iz) * SWICUMCC(j,i-1,k,iz) *
     &            ( ZETA(j,i-1,k,iz) - ZETA(j,i,k,iz) ) +
     &            EPS(iz) * SWICC(j,i-1,k,iz) *
     &            ( ( ZETA(j,i-1,k,iz) - ZETA(j,i-1,k,iz+1) ) -
     &              ( ZETA(j,i,k,iz)   - ZETA(j,i,k,iz+1)   ) )
                RHS(j,i,k) = RHS(j,i,k) - rhoc * switfact
              END IF
C---------------FRONT FACE
              IF (i.NE.NROW) THEN
                rhoc = DELNUS(iz) * SWICUMCC(j,i,k,iz) *
     &            ( ZETA(j,i+1,k,iz) - ZETA(j,i,k,iz) ) +
     &            EPS(iz) * SWICC(j,i,k,iz) *
     &            ( ( ZETA(j,i+1,k,iz) - ZETA(j,i+1,k,iz+1) ) -
     &              ( ZETA(j,i,k,iz)   - ZETA(j,i,k,iz+1)   ) )
                RHS(j,i,k) = RHS(j,i,k) - rhoc * switfact
              END IF
              !LATER=RHS(j,i,k)
!             IF (EERST.NE.LATER) THEN
!               TTTT=3
!             END IF
         END DO IZDDIFF
        END DO JDDIFF
       END DO IDDIFF
      END DO KDDIFF
C
C---------UPDATE ZETA SURFACE BASED ON CURRENT HEAD SOLUTION
C         PRIOR TO SOLVING FOR NEW HEAD SOLUTION IN SELECTED
C         GWF PROCESS SOLVER
        CALL SSWIUPDZ(Kkstp,Kkper,Kkiter)
C
C---------EVALUATE IF THE NUMBER OF SWI TIME STEPS PER MODFLOW
C         TIME STEP SHOULD BE INCREASED AND THE SWI SOLUTION
C         RESTARTED AT THE FIRST SWI TIME STEP IN THE CURRENT
C         MODFLOW TIME STEP.
2010    FORMAT(1X,'IADPT      :',1X,I10,
     2         1X,'ADJ. FACTOR:',1X,G10.3)
        IF ( NADPTFLG.NE.0 ) THEN
          r = REAL( IADPT, 4 )
          IF ( IADPTMOD.GT.0 .AND. IADPT.NE.NADPTMX ) THEN
            WRITE (IOUT,2010) IADPT,ADPTVAL
!            f = ADPTVAL / ADPTFCT
            r = r * ADPTVAL
            IF ( r.GT.REAL( NADPTMX, 4 ) ) THEN
              r = REAL( NADPTMX, 4 )
            END IF
C-------------RESET RHS AND ZETA IF CHANGE IN ZETA EXCEEDS
            DO k = 1, NLAY
              DO i = 1, NROW
                DO j = 1, NCOL
                  RHS(j,i,k)        = RHSFRESH(j,i,k)
                  QZEXTRACUM(j,i,k) = 0.0
                  DO iz = 2, NZONES
!                    ZETA(j,i,k,iz) = ZETAIT(j,i,k,iz)
                    ZETA(j,i,k,iz) = ZETAOLD(j,i,k,iz)
                  END DO
                END DO
              END DO
            END DO
C-------------UPDATE IADPT, IADPTMOD, AND SWIDELT
            IADPT    = CEILING( r )
            IADPTMOD = 0
            switfact = 1.0 / REAL( IADPT, 4 )
            SWIDELT  = DELT * switfact
            WRITE (IOUT,2000) IADPT, DELT,SWIDELT
C------------RESTART SWI CALCULATIONS AT THE BEGINNING
           iswi = 1
           ADPTVAL = 1.0
           CYCLE ADAPTIVE
          END IF
C           RESET IADPTMOD AND ADPTVAL
          IADPTMOD = -1
          ADPTVAL  = 1.0
        END IF
C
C---------MOVE TIPS AND TOES OF ZETA SURFACES
        CALL SSWITIPTOE()
!C
!C---------SAVE ZETA TO ZETAIT IF THIS IS THE FIRST ADAPTIVE
!        IF ( iswi.EQ.1 ) THEN
!          DO k = 1, NLAY
!            DO i = 1, NROW
!              DO j = 1, NCOL
!                DO iz = 2, NZONES
!                  ZETAIT(j,i,k,iz) = ZETA(j,i,k,iz)
!                END DO
!              END DO
!            END DO
!          END DO
!        END IF
C
C-------END OF ADAPTIVE
        iswi = iswi + 1
        IF ( iswi.GT.IADPT ) EXIT ADAPTIVE
      END DO ADAPTIVE
C
C---------RETURN
        RETURN
      END SUBROUTINE GWF2SWI1FM


      SUBROUTINE GWF2SWI1BD(Kkstp,Kkper,Kkiter,Igrid)
C-----VERSION 3 26SEP2010 GWF2SWI1BD
C     ******************************************************************
C     CALCULATE AND SAVE ZETA FOR SWI PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                        LBOTM,BOTM,
     3                        CR,CC,CV,HCOF,RHS,
     4                        DELR,DELC,IBOUND,HNEW,HOLD,
     5                        BUFF,ISSFLG,NSTP
        USE GWFBASMODULE, ONLY: DELT,HDRY,TOTIM,
     2                          IHDDFL,IBUDFL,ICBCFL,IOFLG,
     3                          MSUM,VBVL,VBNM
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Kkstp
        INTEGER, INTENT(IN) :: Kkper
        INTEGER, INTENT(IN) :: Kkiter
        INTEGER, INTENT(IN) :: Igrid
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: iozeta
        INTEGER :: iu
        INTEGER :: ibd, ibdlbl
        INTEGER :: i, j, k, n
        INTEGER :: iz
        REAL :: zero, rate, rin, rout
        DOUBLE PRECISION :: dzero, ratin, ratout, rrate
        REAL :: z
        CHARACTER*16               :: ZETANAME
        CHARACTER*16, DIMENSION(1) :: text
        CHARACTER*15               :: cobs
C
        DATA text /'     SWI QZEXTRA'/
C
C-------OUTPUT FORMAT STATEMENTS
 2000  FORMAT(1X,/1X,A,'   PERIOD ',I4,'   STEP ',I3)
 2010  FORMAT(1X,'   LAYER ',I3,'   ROW ',I5,'   COL ',
     1        I5,'   RATE ',1PG15.6)
 2020  FORMAT(1X,'SOLVING FOR ZETA FOR LAYER',1X,I5,1X,'ZONE',1X,I5)
 2030  FORMAT('    ZETAPLANE ',I2)
 2040  FORMAT(A15,1X)
 2050  FORMAT(A12,I3.3,1X)
 2060  FORMAT(1X)
 2070  FORMAT(G15.7)
C
C-------SET POINTERS FOR GRID
      CALL SGWF2SWI1PNT(Igrid)
C
C-------INITIALIZE ZERO VARIABLES
      zero   = 0.0
      dzero  = 0.0D0
C
C-------DETERMINE IF ZETA VALUES SHOULD BE OUTPUT
      iozeta = 0
      IF ( IHDDFL.GT.0 ) THEN
        DO k = 1, NLAY
!          IF ( IOFLG(k,3).GT.0 ) THEN
          IF ( IOFLG(k,3).NE.0 ) THEN
            iozeta = 1
            EXIT
          END IF
        END DO
      END IF
C
C-------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND
      ibd    = 0
      IF( ISWICB.LT.0 .AND. ICBCFL.NE.0 ) ibd = -1
      IF( ISWICB.GT.0 ) ibd = ICBCFL
      ibdlbl = 0
C
C-------QZEXTRA FLOW TERMS
      ibdlbl = 0
      ratout = dzero
      ratin  = dzero
      DO k=1,NLAY
        DO i=1,NROW
          DO j=1,NCOL
            BUFF(j,i,k) = zero
            IF ( IBOUND(j,i,k).LT.1 ) CYCLE
            rrate = QZEXTRACUM(j,i,k)
            rate  = rrate
C
C-------------PRINT THE INDIVIDUAL RATES IF REQUESTED(IGHBCB<0).
            IF( ibd.LT.0 ) THEN
               IF( ibdlbl.EQ.0 ) WRITE(IOUT,2000) text(1),Kkper,Kkstp
               WRITE(IOUT,2010) k, i, j, rate
               ibdlbl=1
            END IF
C
C-------------ADD RATE TO BUFFER.
            BUFF(j,i,k) = BUFF(j,i,k) + rate
C
C-------------SEE IF FLOW IS INTO AQUIFER OR OUT OF AQUIFER.
            IF( rate.LT.zero ) THEN
C
C---------------FLOW IS OUT OF AQUIFER SUBTRACT RATE FROM RATOUT.
              ratout = ratout - rrate
            ELSE
C
C---------------FLOW IS INTO AQIFER; ADD RATE TO RATIN.
              ratin = ratin + rrate
            END IF
          END DO
        END DO
      END DO
C
C-------IF CELL-BY-CELL TERMS WILL BE SAVED AS A 3-D ARRAY, THEN CALL
C       UTILITY MODULE UBUDSV TO SAVE THEM.
      IF ( ibd.EQ.1 ) CALL UBUDSV(Kkstp,Kkper,text(1),ISWICB,BUFF,
     2                            NCOL,NROW,NLAY,IOUT)
!C
!C-------MOVE RATES, VOLUMES AND LABELS INTO ARRAYS FOR PRINTING.
!      rin          = ratin
!      rout         = ratout
!      VBVL(3,MSUM) = rin
!      VBVL(1,MSUM) = VBVL(1,MSUM) + rin  * DELT
!      VBVL(4,MSUM) = rout
!      VBVL(2,MSUM) = VBVL(2,MSUM) + rout * DELT
!      VBNM(MSUM)   = text(1)
!C
!C-------INCREMENT THE BUDGET TERM COUNTER.
!      MSUM = MSUM + 1
!C
!C---------UPDATE ZETA SURFACE BASED ON FINAL HEAD SOLUTION 
!C         FOR MODFLOW TIMESTEP
!        CALL SSWIUPDZ(Kkstp,Kkper,Kkiter)
!C
!C---------MOVE TIPS AND TOES OF ZETA SURFACES
!        CALL SSWITIPTOE()
C
C---------WRITE ZETA TO UNFORMATTED FILE
!        IF ( iozeta.GT.0 .AND. ISWIZT.NE.0 ) THEN
        IF ( iozeta.NE.0 .AND. ISWIZT.NE.0 ) THEN
          DO iz = 2, NZONES
            WRITE(ZETANAME,2030) iz-1
            CALL UBUDSV(kkstp,Kkper,ZETANAME,
     &                  ISWIZT,ZETA(1,1,1,iz),NCOL,NROW,NLAY,IOUT)
          END DO
        END IF
C
C---------WRITE OBSERVATION DATA TO OBSERVATION FILE
!        IF ( NOBS.GT.0 .AND. ISWIOBS.NE.0 ) THEN
        IF ( NOBS.NE.0 .AND. ISWIOBS.NE.0 ) THEN
          iu = ABS( ISWIOBS )
C           WRITE HEADER FOR ASCII FILE
          IF ( ISWIOBS.GT.0 ) THEN
            IF ( IOBSHEADER.LT.1 ) THEN
              WRITE (iu,2040,ADVANCE='NO') '          TOTIM'
              DO n = 1, NOBS
                DO iz = 2, NZONES
                  WRITE (iu,2050,ADVANCE='NO')
     2              ADJUSTR(TRIM(SWIOBS(n)%OBSNAM)), iz-1
                END DO
              END DO
              WRITE (iu,2060,ADVANCE='YES')
              IOBSHEADER = 1
            END IF
          END IF
C           WRITE DATA TO OBSERVATION OUTPUT FILE
          BINOBS: DO n = 1, NOBS
            k = SWIOBS(n)%KLAY
            i = SWIOBS(n)%IROW
            j = SWIOBS(n)%JCOL
            IZBINOBS: DO iz = 2, NZONES
              z = ZETA(j,i,k,iz)
C               ASCII OUTPUT
              IF ( ISWIOBS.GT.0 ) THEN
                IF ( n.EQ.1 .AND. iz.EQ.2 ) THEN
                  WRITE (cobs,2070) TOTIM
                  WRITE (iu,2040,ADVANCE='NO') ADJUSTR(TRIM(cobs))
                END IF
                WRITE (cobs,2070) z
                WRITE (iu,2040,ADVANCE='NO') ADJUSTR(TRIM(cobs))
                IF ( n.EQ.NOBS .AND. iz.EQ.NZONES) THEN
                  WRITE (iu,2060,ADVANCE='YES')
                END IF
C               BINARY OUTPUT
              ELSE
                IF ( n.EQ.1 .AND. iz.EQ.2 ) WRITE (iu) TOTIM
                WRITE (iu) z
              END IF
            END DO IZBINOBS
          END DO BINOBS
        END IF
C
C---------RETURN
        RETURN
      END SUBROUTINE GWF2SWI1BD


      SUBROUTINE SSWI1SQR(J,I,K,IZ,NCOL,NROW,NLAY,NZONES,
     &                   Q,HNEW,ZETA,DELNUS,EPS,SWICUMCR,SWICR)

C     ******************************************************************
C     SWI - CALCULATE FLUX IN THE ROW DIRECTION
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: J,I,K,IZ
        INTEGER, INTENT(IN) :: NCOL,NROW,NLAY,NZONES
        REAL, INTENT(INOUT) :: Q
        DOUBLEPRECISION, DIMENSION(NCOL,NROW,NLAY), INTENT(IN) :: HNEW
        REAL, DIMENSION(NCOL,NROW,NLAY,NZONES+1), INTENT(IN) :: ZETA
        REAL, DIMENSION(NZONES), INTENT(IN) :: DELNUS
        REAL, DIMENSION(NZONES), INTENT(IN) :: EPS
        REAL, DIMENSION(NCOL,NROW,NLAY,NZONES), INTENT(IN)   :: SWICR
        REAL, DIMENSION(NCOL,NROW,NLAY,NZONES), INTENT(IN)   :: SWICUMCR
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: iz2
C
C       + + + CODE + + +
C
C---------CALCULATE Q IN ROW DIRECTION
        Q = SWICUMCR(J-1,I,K,IZ) * (HNEW(J-1,I,K)-HNEW(J,I,K))

        DO IZ2 = 1, IZ
          Q = Q + SWICUMCR(J-1,I,K,IZ) * DELNUS(IZ) *
     &            (ZETA(J-1,I,K,IZ2)-ZETA(J,I,K,IZ2))
        END DO

        Q = Q - SWICR(J-1,I,K,IZ) * EPS(IZ) *
     &          (ZETA(J-1,I,K,IZ)-ZETA(J,I,K,IZ))

        Q = Q + SWICR(J-1,I,K,IZ) * EPS(IZ) *
     &          (ZETA(J-1,I,K,IZ+1)-ZETA(J,I,K,IZ+1))


        DO IZ2 = IZ+1, NZONES
          Q = Q + SWICUMCR(J-1,I,K,IZ2) * DELNUS(IZ2) *
     &            (ZETA(J-1,I,K,IZ2)-ZETA(J,I,K,IZ2))
          Q = Q - SWICR(J-1,I,K,IZ2) * EPS(IZ2) *
     &            (ZETA(J-1,I,K,IZ2)-ZETA(J,I,K,IZ2))
          Q = Q + SWICR(J-1,I,K,IZ2) * EPS(IZ2) *
     &            (ZETA(J-1,I,K,IZ2+1)-ZETA(J,I,K,IZ2+1))
        END DO
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWI1SQR



      SUBROUTINE SSWI1SQC(J,I,K,IZ,NCOL,NROW,NLAY,NZONES,
     &                   Q,HNEW,ZETA,DELNUS,EPS,SWICUMCC,SWICC)
C
C     ******************************************************************
C     SWI - CALCULATE FLUX IN THE COLUMN DIRECTION
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: J,I,K,IZ
        INTEGER, INTENT(IN) :: NCOL,NROW,NLAY,NZONES
        REAL, INTENT(INOUT) :: Q
        DOUBLEPRECISION, DIMENSION(NCOL,NROW,NLAY), INTENT(IN) :: HNEW
        REAL, DIMENSION(NCOL,NROW,NLAY,NZONES+1), INTENT(IN) :: ZETA
        REAL, DIMENSION(NZONES), INTENT(IN) :: DELNUS
        REAL, DIMENSION(NZONES), INTENT(IN) :: EPS
        REAL, DIMENSION(NCOL,NROW,NLAY,NZONES), INTENT(IN)   :: SWICC
        REAL, DIMENSION(NCOL,NROW,NLAY,NZONES), INTENT(IN)   :: SWICUMCC
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: iz2
C
C       + + + CODE + + +
C
C---------CALCULATE Q IN COLUMN DIRECTION
        Q = SWICUMCC(J,I-1,K,IZ) * (HNEW(J,I-1,K)-HNEW(J,I,K))
C
        DO IZ2 = 1, IZ
          Q = Q + SWICUMCC(J,I-1,K,IZ) * DELNUS(IZ) *
     &            (ZETA(J,I-1,K,IZ2)-ZETA(J,I,K,IZ2))
        END DO
C
        Q = Q - SWICC(J,I-1,K,IZ) * EPS(IZ) *
     &          (ZETA(J,I-1,K,IZ)-ZETA(J,I,K,IZ))
        Q = Q + SWICC(J,I-1,K,IZ) * EPS(IZ) *
     &          (ZETA(J,I-1,K,IZ+1)-ZETA(J,I,K,IZ+1))

C
        DO IZ2 = IZ + 1, NZONES
          Q = Q + SWICUMCC(J,I-1,K,IZ2) * DELNUS(IZ2) *
     &            (ZETA(J,I-1,K,IZ2)-ZETA(J,I,K,IZ2))
          Q = Q - SWICC(J,I-1,K,IZ2) * EPS(IZ2) *
     &            (ZETA(J,I-1,K,IZ2)-ZETA(J,I,K,IZ2))
          Q = Q + SWICC(J,I-1,K,IZ2) * EPS(IZ2) *
     &            (ZETA(J,I-1,K,IZ2+1)-ZETA(J,I,K,IZ2+1))
        END DO
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWI1SQC


      SUBROUTINE SSWI1SD(J,I,NCOL,NROW,IPOS,FAC,B,CRLAY,CCLAY,VAR)
C
C     ******************************************************************
C     SWI - CALCULATE FLUX IN THE COLUMN AND ROW DIRECTIONS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: J,I
        INTEGER, INTENT(IN) :: NCOL,NROW
        INTEGER, DIMENSION(NCOL,NROW), INTENT(IN) :: IPOS
        REAL, INTENT(IN) :: FAC
        REAL, DIMENSION(NCOL,NROW), INTENT(INOUT) :: B
        REAL, DIMENSION(NCOL,NROW), INTENT(IN) :: CRLAY
        REAL, DIMENSION(NCOL,NROW), INTENT(IN) :: CCLAY
        DOUBLEPRECISION, DIMENSION(NCOL,NROW), INTENT(IN) :: VAR
C       + + + LOCAL DEFINITIONS + + +
C       + + + CODE + + +
C
C---------CALCULATE FLUX IN COLUMN DIRECTION
C         (LEFT FACE)
        IF ( J.NE.1 ) THEN
          IF (IPOS(J-1,I).EQ.0) THEN
            B(J,I) = B(J,I)- FAC*CRLAY(J-1,I)*(VAR(J-1,I)-VAR(J,I))
          END IF
        END IF
C         (RIGHT FACE)
        IF ( J.NE.NCOL ) THEN
          IF (IPOS(J+1,I).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CRLAY(J,I)*(VAR(J+1,I)-VAR(J,I))
          END IF
        END IF
C
C---------CALCULATE FLUX IN ROW DIRECTION
C         (BACK FACE)
        IF ( I.NE.1 ) THEN
          IF (IPOS(J,I-1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I-1)*(VAR(J,I-1)-VAR(J,I))
          END IF
        END IF
C         (FRONT FACE)
        IF ( I.NE.NROW ) THEN
          IF (IPOS(J,I+1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I)*(VAR(J,I+1)-VAR(J,I))
          END IF
        END IF
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWI1SD


      SUBROUTINE SSWI1SR(J,I,NCOL,NROW,IPOS,FAC,B,CRLAY,CCLAY,VAR)
C
C     ******************************************************************
C     SWI - CALCULATE FLUX IN THE COLUMN AND ROW DIRECTIONS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: J,I
        INTEGER, INTENT(IN) :: NCOL,NROW
        INTEGER, DIMENSION(NCOL,NROW), INTENT(IN) :: IPOS
        REAL, INTENT(IN) :: FAC
        REAL, DIMENSION(NCOL,NROW), INTENT(INOUT) :: B
        REAL, DIMENSION(NCOL,NROW), INTENT(IN) :: CRLAY
        REAL, DIMENSION(NCOL,NROW), INTENT(IN) :: CCLAY
        REAL, DIMENSION(NCOL,NROW), INTENT(IN) :: VAR
C       + + + LOCAL DEFINITIONS + + +
C       + + + CODE + + +
C
C----------CALCULATE FLUX IN COLUMN DIRECTION
C         (LEFT FACE)
        IF ( J.NE.1 ) THEN
          IF (IPOS(J-1,I).EQ.0) THEN
            B(J,I) = B(J,I)- FAC*CRLAY(J-1,I)*(VAR(J-1,I)-VAR(J,I))
          END IF
        END IF
C         (RIGHT FACE)
        IF ( J.NE.NCOL ) THEN
          IF (IPOS(J+1,I).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CRLAY(J,I)*(VAR(J+1,I)-VAR(J,I))
          END IF
        END IF

C
C--------CALCULATE FLUX IN ROW DIRECTION
C        (BACK FACE)
        IF ( I.NE.1 ) THEN
          IF (IPOS(J,I-1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I-1)*(VAR(J,I-1)-VAR(J,I))
          END IF
        END IF
C        (FRONT FACE)
        IF ( I.NE.NROW ) THEN
          IF (IPOS(J,I+1).EQ.0) THEN
            B(J,I) = B(J,I)-FAC*CCLAY(J,I)*(VAR(J,I+1)-VAR(J,I))
          END IF
        END IF
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWI1SR

      SUBROUTINE SSWIUPDZ(Kkstp,Kkper,Kkiter)
C-----VERSION 3 26SEP2010 SSWIUPDZ
C     ******************************************************************
C     UPDATE ZETA SURFACES FOR SWI PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                        LBOTM,BOTM,
     3                        CR,CC,CV,HCOF,RHS,
     4                        DELR,DELC,IBOUND,HNEW,HOLD,
     5                        BUFF,ISSFLG,NSTP
        USE GWFBASMODULE, ONLY: DELT,HDRY,TOTIM
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Kkstp
        INTEGER, INTENT(IN) :: Kkper
        INTEGER, INTENT(IN) :: Kkiter
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k, n
        INTEGER :: kk
        INTEGER :: iad
        INTEGER :: iz, iz2, iz3
        INTEGER :: izrev
        INTEGER :: iusezone, iactive
        INTEGER :: icount
        INTEGER :: mxiterde45, niter
        INTEGER :: ksoliter, kksoliter
        INTEGER :: iprd4, ifreq, ierr, icnvg
        REAL :: z, z0, dz, f, r
        REAL :: q, qztop, qzbot
        REAL :: accl, hclose
        REAL :: one, tled, rho
        REAL :: strg
        REAL :: zt, zb, zetac, zcalc, zdiff, zetaavg
        REAL :: nuontop, nubelbot
        REAL :: zclose, rclose
C
C-------OUTPUT FORMAT STATEMENTS
 2000  FORMAT(1X,/1X,A,'   PERIOD ',I4,'   STEP ',I3)
 2010  FORMAT(1X,'   LAYER ',I3,'   ROW ',I5,'   COL ',
     1        I5,'   RATE ',1PG15.6)
 2020  FORMAT(1X,'SOLVING FOR ZETA FOR LAYER',1X,I5,1X,'ZONE',1X,I5)
 2030  FORMAT('    ZETAPLANE ',I2)
C
C---------CALCULATE BRHS
        KBRHS: DO k=1,NLAY
          IBRHS: DO i=1,NROW
            JBRHS: DO j=1,NCOL
              IZBRHS: DO iz=1,NZONES
C
C                 INITIALIZE BRHS
                BRHS(j,i,iz)=0.
C
C                 CALCULATE BOUNDARY FLUX
                q = RHSFRESH(j,i,k) -
     2              HNEW(j,i,k)*(HCOF(j,i,k)-HCOFFRESH(j,i,k))

                IF (ISSFLG(Kkper).NE.0)  THEN
                  strg = 0.
!                  iusezone=IZONENR(j,i,k)
!                   IF ((IZONENR(j,i,k).LT.0).AND.(q.GT.0)) THEN
!                           iusezone=1
!                   ENDIF
                ELSE
                  one  = 1.
!                  tled = one/DELT
                  tled = one / SWIDELT
                  rho  = SC1(j,i,k) * tled
                  q    = q + rho * HOLD(j,i,k)
!                  iusezone=IZONENR(j,i,k)
                  strg = rho*HOLD(j,i,k) - rho*HNEW(j,i,k)
                  q = q - strg
!                  IF ((IZONENR(j,i,k).LT.0).AND.(q.GT.0)) THEN
!                    iusezone=1
!                  ENDIF
                END IF
C
C                 SET ZONE NUMBER FOR BOUNDARY CONDITIONS
                iusezone=IZONENR(j,i,k)
                IF ( (IZONENR(j,i,k).LT.0) .AND. (q.GT.0) ) THEN
                  iusezone = 1
                ENDIF

                IF (IPLPOS(j,i,k,iz).EQ.0) THEN
                  IF ((iz.LE.ABS(iusezone)).AND.(q.NE.0)) THEN
                    BRHS(j,i,iz)=BRHS(j,i,iz) + q
                  ENDIF
                ELSE
                  BRHS(j,i,iz) = 0.0
                ENDIF
C
                IF ( iz.NE.1 ) THEN
                  IF (IPLPOS(j,i,k,iz).EQ.0) THEN
!                    BRHS(j,i,iz)= BRHS(j,i,iz) -
!     &                SSZ(j,i,k)*DELR(j)*DELC(i)*ZETA(j,i,k,iz)/DELT
                    BRHS(j,i,iz)= BRHS(j,i,iz) -
     &                SSZ(j,i,k)*DELR(j)*DELC(i)*ZETA(j,i,k,iz)/SWIDELT
                    CALL SSWI1SD(j,i,NCOL,NROW,IPLPOS(1,1,k,iz),
     &                1.,BRHS(1,1,iz),
     &                SWICUMCR(1,1,k,iz),SWICUMCC(1,1,k,iz),
     &                HNEW(1,1,k))
                    DO iz2 = 1, iz-1
                      CALL SSWI1SR(j,i,NCOL,NROW,IPLPOS(1,1,k,iz),
     &                  DELNUS(iz2),BRHS(1,1,iz),
     &                  SWICUMCR(1,1,k,iz),SWICUMCC(1,1,k,iz),
     &                  ZETA(1,1,k,iz2))
                    END DO
                  END IF
                  CALL SSWI1SR(j,i,NCOL,NROW,IPLPOS(1,1,k,iz),
     &              EPS(iz),BRHS(1,1,iz),
     &              SWICR(1,1,k,iz),SWICC(1,1,k,iz),
     &              ZETA(1,1,k,iz+1))
C
C-----------Zones iz+1 through bottom
                  DO iz2 = iz+1, NZONES
                    IF (IPLPOS(j,i,k,iz).EQ.0) THEN
                      CALL SSWI1SR(j,i,NCOL,NROW,IPLPOS(1,1,k,iz),
     &                  DELNUS(iz2),BRHS(1,1,iz),
     &                  SWICUMCR(1,1,k,iz2),SWICUMCC(1,1,k,iz2),
     &                  ZETA(1,1,k,iz2))
                    END IF
C---------------------LEFT FACE
                    IF ( j.NE.1 ) THEN
                      BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &                  SWICR(j-1,i,k,iz2)*
     &                  ((ZETA(j-1,i,k,iz2)-ZETA(j,i,k,iz2))-
     &                  (ZETA(j-1,i,k,iz2+1)-ZETA(j,i,k,iz2+1)))
                    END IF
C---------------------RIGHT FACE
                    IF ( j.NE.NCOL ) THEN
                      BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &                  SWICR(j,i,k,iz2)*
     &                 ((ZETA(j+1,i,k,iz2)-ZETA(j,i,k,iz2))-
     &                 (ZETA(j+1,i,k,iz2+1)-ZETA(j,i,k,iz2+1)))
                    END IF
C---------------------BACK FACE
                    IF ( i.NE.1 ) THEN
                      BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &                  SWICC(j,i-1,k,iz2)*
     &                  ((ZETA(j,i-1,k,iz2)-ZETA(j,i,k,iz2))-
     &                  (ZETA(j,i-1,k,iz2+1)-ZETA(j,i,k,iz2+1)))
                    END IF
C---------------------FRONT FACE
                    IF ( i.NE.NROW ) THEN
                      BRHS(j,i,iz) = BRHS(j,i,iz)+EPS(iz2)*
     &                  SWICC(j,i,k,iz2)*
     &                  ((ZETA(j,i+1,k,iz2)-ZETA(j,i,k,iz2))-
     &                  (ZETA(j,i+1,k,iz2+1)-ZETA(j,i,k,iz2+1)))
                    END IF
                  END DO
C
C-------------------SPECIFY BOUNDARY CONDITION OF FLOW AT TIPS AND TOES
                  IF (IPLPOS(j,i,k,iz).EQ.0) THEN
C---------------------LEFT FACE
                    IF ( j.NE.1 ) THEN
                      IF (IPLPOS(j-1,i,k,iz).NE.0) THEN
                        CALL SSWI1SQR(j,i,k,iz,NCOL,NROW,NLAY,NZONES,
     &                         q,HNEW,ZETA,DELNUS,EPS,
     &                         SWICUMCR,SWICR)
                        BRHS(j,i,iz) = BRHS(j,i,iz) - q
                      END IF
                    END IF
C---------------------RIGHT FACE
                    IF ( j.NE.NCOL ) THEN
                      IF (IPLPOS(j+1,i,k,iz).NE.0) THEN
                        CALL SSWI1SQR(j+1,i,k,iz,NCOL,NROW,NLAY,NZONES,
     &                         q,HNEW,ZETA,DELNUS,EPS,
     &                         SWICUMCR,SWICR)
                        BRHS(j,i,iz) = BRHS(j,i,iz) + q
                      END IF
                    END IF
C---------------------BACK FACE
                    IF ( i.NE.1 ) THEN
                      IF (IPLPOS(j,i-1,k,iz).NE.0) THEN
                        CALL SSWI1SQC(j,i,k,iz,NCOL,NROW,NLAY,NZONES,
     &                         q,HNEW,ZETA,DELNUS,EPS,
     &                         SWICUMCC,SWICC)
                        BRHS(j,i,iz) = BRHS(j,i,iz) - q
                      END IF
                    END IF
C---------------------FRONT FACE
                    IF ( i.NE.NROW ) THEN
                      IF (IPLPOS(j,i+1,k,iz).NE.0) THEN
                        CALL SSWI1SQC(j,i+1,k,iz,NCOL,NROW,NLAY,NZONES,
     &                         q,HNEW,ZETA,DELNUS,EPS,
     &                         SWICUMCC,SWICC)
                        BRHS(j,i,iz) = BRHS(j,i,iz) + q
                      END IF
                    END IF
C---------------------UPPER FACE
                    IF ( k.NE.1 ) THEN
                      qztop = CV(j,i,k-1)*(HNEW(j,i,k)-HNEW(j,i,k-1))+
     &                        QZEXTRA(j,i,k)
                      nuontop = NUBOT(j,i,k-1)
                      IF ((qztop.LT.0).AND.(nuontop.GE.NUS(iz)).AND.
     &                    (NUBOT(j,i,k).GE.nuontop)) THEN
                        BRHS(j,i,iz) = BRHS(j,i,iz) + qztop
                      ENDIF
                    END IF
C---------------------LOWER FACE
                    IF ( k.NE.NLAY ) THEN
                      qzbot = CV(j,i,k)*(HNEW(j,i,k)-HNEW(j,i,k+1))-
     &                        QZEXTRA(j,i,k+1)
                      nubelbot = NUTOP(j,i,k+1)
                      IF ((qzbot.LT.0).AND.(nubelbot.LT.NUS(iz)).AND.
     &                    (NUTOP(j,i,k).LE.nubelbot)) THEN
                        BRHS(j,i,iz) = BRHS(j,i,iz) + 0
                      ELSE
                        BRHS(j,i,iz) = BRHS(j,i,iz) + qzbot
                      END IF
                    END IF
                  END IF
                END IF

C               BUFF(j,i,iz)=BRHS(j,i,iz)

              END DO IZBRHS
            END DO JBRHS
          END DO IBRHS
C
C-----------SOLVE EQUATIONS FOR 1 LAYER
          IZSOLVZ4LAY: DO iz = 2, NZONES
            iactive = 0
C             CALCULATE SWI CONDUCTANCES FOR CURRENT ZETA SURFACE
C             IN CURRENT LAYER
            ISOLVZ4LAY: DO i=1,NROW
              JSOLVZ4LAY: DO j=1,NCOL
C
C                 CALCULATE SWISOLVCR
                IF ( j.NE.NCOL ) THEN
                  IF ((IPLPOS(j,i,k,iz).EQ.0).AND.
     &                (IPLPOS(j+1,i,k,iz).EQ.0)) THEN
                    SWISOLVCR(j,i) = DELNUS(iz)*SWICUMCR(j,i,k,iz)
                  ELSE
                    SWISOLVCR(j,i) = 0.
                  END IF
                END IF
                SWISOLVCR(j,i) = SWISOLVCR(j,i)-EPS(iz)*SWICR(j,i,k,iz)
C
C                 CALCULATE SWISOLVCC
                IF ( i.NE.NROW ) THEN
                  IF ((IPLPOS(j,i,k,iz).EQ.0).AND.
     &                (IPLPOS(j,i+1,k,iz).EQ.0)) THEN
                    SWISOLVCC(j,i) = DELNUS(iz)*SWICUMCC(j,i,k,iz)
                  ELSE
                    SWISOLVCC(j,i)=0.
                  END IF
                END IF
                SWISOLVCC(j,i) = SWISOLVCC(j,i)-EPS(iz)*SWICC(j,i,k,iz)
C
C                 CALCULATE SWISOLVCV AND SWIHCOF (STORAGE)
                SWISOLVCV(j,i) = 0.
!                SWIHCOF(j,i)   = -SSZ(j,i,k)*DELR(j)*DELC(i)/DELT
                SWIHCOF(j,i)   = -SSZ(j,i,k)*DELR(j)*DELC(i)/SWIDELT

                DUM(j,i)       = ZETA(j,i,k,iz)
C
C                 DETERMINE IF CURRENT CELL HAS ZETA VALUE GREATER THAN
C                 THE BOTTOM AND LESS THAN THE TOP (IPLPOS=0)
                IF (IPLPOS(j,i,k,iz).EQ.0) THEN
                  IBO(j,i) = 1
                  iactive = iactive + 1
                ELSE
                  IBO(j,i) = 0
                END IF
              END DO JSOLVZ4LAY
            END DO ISOLVZ4LAY
C
C-------------SOLVE FOR ZETA SURFACE
            IF ( iactive.NE.0 ) WRITE (IOUT,2020) k, iz
C             DIRECT SOLVER
            SOLVEZETA: IF (NSOLVER.EQ.1) THEN
C-------------INITIALIZE PARAMETERS FOR DIRECT SOLVING
              zclose    = 0.0
              kksoliter = 1
              ierr      = 0
C---------------CALL DE45AP TO SOLVE FOR ZETA
              ACTIVEDE4: IF ( iactive.NE.0 ) THEN
                !IACTIVEEQ1: IF ( iactive.EQ.1 ) THEN
                  DO i=1,NROW
                    DO j=1,NCOL
                      IF ( IBO(j,i).EQ.1 ) THEN
                        qzbot = BRHS(j,i,iz) / SWIHCOF(j,i)
                        !ZETA(j,i,k,iz) = qzbot
                      END IF
                    END DO
                  END DO
                !ELSE
                  CALL DE47AP(DUM,IBO,SWIDE4%AU,SWIDE4%AL,SWIDE4%IUPPNT,
     2              SWIDE4%IEQPNT,SWIDE4%D4B,SWIDE4%MXUP,SWIDE4%MXLOW,
     3              SWIDE4%MXEQ,SWIDE4%MXBW,
     4              SWISOLVCR,SWISOLVCC,SWISOLVCV,SWIHCOF,
     5              BRHS(1,1,iz),SWIDE4%ACCLDE4,kksoliter,
     6              SWIDE4%ITMX,SWIDE4%MXITER,SWIDE4%NITERDE4,
     7              zclose,IPRSOL,icnvg,NCOL,NROW,NLAYSWI,
     8              IOUT,SWIDE4%LRCHDE4,SWIDE4%HDCGDE4,SWIDE4%IFREQ,
     9              Kkstp,Kkper,SWIDE4%DELTL,NSTP(Kkper),SWIDE4%ID4DIR,
     x              SWIDE4%ID4DIM,MUTSOL,SWIDE4%DELTL,SWIDE4%NBWL,
     1              SWIDE4%NUPL,SWIDE4%NLOWL,SWIDE4%NLOW,SWIDE4%NEQ,
     2              SWIDE4%NUP,SWIDE4%NBW,ierr)

                  IF (ierr.EQ.0) THEN
                    WRITE(IOUT,'(A,I3,A,I3,A,I3,A,I3,A,I5,A,I8)')
     &               ' STRESS PERIOD=',Kkper,
     &               ' TIMESTEP=',kkstp,
     &               ' LAYER=',k,
     &               ' ZONE=',iz,
     &               ' ACTIVE CELLS=',iactive,
     &               ' NUMBER OF EQUATIONS=',SWIDE4%NEQ  !ierr
                  END IF
                !END IF IACTIVEEQ1
              END IF ACTIVEDE4
C             PRECONDITIONED CONJUGATE GRADIENT
            ELSE IF ( NSOLVER.EQ.2 ) THEN
              kksoliter=1
              ierr = 0
C---------------CALL PCG7AP TO SOLVE FOR ZETA
              ACTIVEPCG: IF ( iactive.NE.0 ) THEN
C
C-----------------REINITIALIZE BUFF - USED FOR RES IN PCG
                DO kk = 1, NLAYSWI
                  DO i = 1, NROW
                    DO j = 1, NCOL
                      BUFF(j,i,kk) = 0.0
                    END DO
                  END DO
                END DO
C
C                 LINEAR EQUATION -- OUTER ITERATIONS ONLY NEEDED IF
C                 HCLOSE AND/OR RCLOSE NOT ACHIEVED AFTER FIRST
C                 INNER ITERATION
                zclose = SWIPCG%ZCLOSEPCG
                rclose = ( 1.0 + DELNUS(iz) ) * SWIPCG%RCLOSEPCG
                OUTERPCG: DO ksoliter = 1, SWIPCG%MXITER
                  kksoliter = ksoliter
                  CALL PCG7AP(DUM,IBO,SWISOLVCR,SWISOLVCC,SWISOLVCV,
     2              SWIHCOF,BRHS(1,1,iz),SWIPCG%VPCG,SWIPCG%SS,SWIPCG%P,
     3              SWIPCG%CD,SWIPCG%HCHG,SWIPCG%LHCH,SWIPCG%RCHG,
     4              SWIPCG%LRCHPCG,kksoliter,SWIPCG%NITER,zclose,
     5              rclose,icnvg,Kkstp,Kkper,IPRSOL,
     6              SWIPCG%MXITER,SWIPCG%ITER1,SWIPCG%NPCOND,
     7              SWIPCG%NBPOL,NSTP(Kkper),NCOL,NROW,NLAYSWI,
     7              SWIPCG%NODES,SWIPCG%RELAXPCG,IOUT,MUTSOL,
     8              SWIPCG%IT1,SWIPCG%DAMPPCG,BUFF,SWIPCG%HCSV,ierr,
     9              SWIPCG%HPCG,SWIPCG%DAMPPCGT,ISSFLG(Kkper),HDRY,
     X              SWIPCG%IHCOFADD)
                  IF ( icnvg.EQ.1 ) EXIT OUTERPCG
                END DO OUTERPCG
	        END IF ACTIVEPCG
C             UNSTRUCTURED PRECONDITIONED CONJUGATE GRADIENT
            ELSE IF ( NSOLVER.EQ.3 ) THEN
              kksoliter=1
              ierr = 0
C---------------CALL PCGU7AP TO SOLVE FOR ZETA
              ACTIVEPCGU: IF ( iactive.NE.0 ) THEN
C                 LINEAR EQUATION -- OUTER ITERATIONS ONLY NEEDED IF
C                 HCLOSE AND/OR RCLOSE NOT ACHIEVED AFTER FIRST
C                 INNER ITERATION
                zclose = SWIPCGU%ZCLOSEPCGU
                rclose = ( 1.0 + DELNUS(iz) ) * SWIPCGU%RCLOSEPCGU
                OUTERPCGU: DO ksoliter = 1, SWIPCGU%MXITER
                  kksoliter = ksoliter
                   CALL PCGU7AP(icnvg,Kkstp,Kkper,NSTP(Kkper),
     2               SWIPCGU%MXITER,kksoliter,
     3               NCOL,NROW,NLAYSWI,SWIPCGU%NODES,HDRY,IOUT,
     4               SWIPCGU%NPC,SWIPCGU%NORD,
     5               SWIPCGU%ITER1C,SWIPCGU%NITERC,
     6               SWIPCGU%NNZC,SWIPCGU%NIAC,SWIPCGU%NIAPC,
     7               SWIPCGU%NIWC,SWIPCGU%NLAPC,SWIPCGU%NUAPC,
     8               SWIPCGU%NNZAPROC,SWIPCGU%NIAPROC,
     9               zclose,rclose,
     X               DUM,IBO,SWISOLVCR,SWISOLVCC,SWISOLVCV,
     1               SWIHCOF,BRHS(1,1,iz),
     2               IPRSOL,MUTSOL,
     3               SWIPCGU%LHCHPCGU,SWIPCGU%LRCHPCGU,SWIPCGU%IT1PCGU,
     4               SWIPCGU%HCHGPCGU,SWIPCGU%RCHGPCGU,
     5               SWIPCGU%NODEC,
     6               SWIPCGU%BC,SWIPCGU%XC,SWIPCGU%AC,
     7               SWIPCGU%TAPC,SWIPCGU%APC,
     8               SWIPCGU%IAC,SWIPCGU%JAC,SWIPCGU%IUC,
     9               SWIPCGU%IXMAP,SWIPCGU%IWC,
     X               SWIPCGU%IAPC,SWIPCGU%JAPC,
     1               SWIPCGU%IRAPC,SWIPCGU%IAPCMAP,
     2               SWIPCGU%DC,SWIPCGU%PC,SWIPCGU%QC,SWIPCGU%ZC,
     3               SWIPCGU%SCL)
                  IF ( icnvg.EQ.1 ) EXIT OUTERPCGU
                END DO OUTERPCGU
	        END IF ACTIVEPCGU

            END IF SOLVEZETA
C
C-------------COPY REAL ARRAY DUM INTO DOUBLE ARRAY ZETA
!            IF (.NOT.((NSOLVER.EQ.1).AND.(IACTIVE.EQ.1))) then
            IF ( iactive.NE.0 ) THEN
              DO i=1,NROW
                DO j=1,NCOL
                  IF ( IBO(j,i).LT.1 ) CYCLE
                  z  = DUM(j,i)
C                   UPDATE ZETA WITH CURRENT SOLUTION
                  ZETA(j,i,k,iz)   = z
                END DO
              END DO
            END IF
          END DO IZSOLVZ4LAY
        END DO KBRHS
!C
!C---------MOVE TIPS AND TOES OF ZETA SURFACES
!        CALL SSWITIPTOE()
C
C
C---------CHECK SLOPE AT TIPS AND TOES OF ZETA SURFACES
        IF ( NADPTFLG.NE.0 ) CALL SSWICHKSLOPE()
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWIUPDZ
C
C---------MOVE TIPS AND TOES
C---------NOTE: MAY DIFFER FROM MATLAB PROTOTYPE BECAUSE ZETA IS
C---------ALTERED DURING THE SEARCH PROCESS
      SUBROUTINE SSWICHKSLOPE()
C-----VERSION 3 26SEP2010 SSWITIPTOE
C     ******************************************************************
C     MOVE TIPS AND TOES FOR SWI PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                        LBOTM,BOTM,
     3                        CR,CC,CV,HCOF,RHS,
     4                        DELR,DELC,IBOUND,HNEW,HOLD,
     5                        BUFF,ISSFLG,NSTP
        USE GWFBASMODULE, ONLY: DELT,HDRY,TOTIM,IHDDFL,IBUDFL
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k, n
        INTEGER :: iz, iz2, iz3
        INTEGER :: izrev
        INTEGER :: icount
        REAL :: z, zt, zb, zetac, zcalc, zdiff, zetaavg, qztop
        REAL :: rvl, dist, b
          IADPTMOD = 0
          rvl      = 1.0 / ADPTFCT
          ADPTVAL  = 1.0  
          ZTIPTOE: DO k = 1, NLAY
            IZTIPTOE: DO iz = 2, NZONES
              ITIPTOE: DO i = 1, NROW
                JTIPTOE: DO j = 1, NCOL
C
C-------------------SET TEMPORARY ZETA VARIABLES
                  zt    = ZETA(j,i,k,1)
                  zb    = ZETA(j,i,k,NZONES+1)
                  zetac = ZETA(j,i,k,iz)
C
C-------------------ONLY EVALUATE TIP AND TOE MOVEMENT TO ADJACENT CELLS FOR CELLS WHERE
C                   ZETA SURFACE IS NOT AT THE TOP (ILPOS=1) OR BOTTOM (ILPOS=2) OF THE LAYER
                  LACTIVE: IF (IPLPOS(j,i,k,iz).EQ.0) THEN
C--------------------CHECK TIPS AND TOES FOR EACH COLUMN BUT THE FIRST AND LAST
                    LCOL: IF ((j.NE.1).AND.(j.NE.NCOL)) THEN
C----------------------LEFT FACE
                      IF (IPLPOS(j-1,i,k,iz).EQ.1) THEN
                        dist  = 0.5 * ( DELR(j) + DELR(j-1) )
                        zcalc = rvl * TIPSLOPE * dist
                        b     = zt-zetac
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      ELSEIF (IPLPOS(j-1,i,k,iz).EQ.2) THEN
                        dist  = 0.5 * ( DELR(j) + DELR(j-1) )
                        zcalc = rvl * TOESLOPE * dist
                        b     = zetac - zb
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      END IF
C----------------------RIGHT FACE
                      IF (IPLPOS(j+1,i,k,iz).EQ.1) THEN
                        dist  = 0.5 * ( DELR(j) + DELR(j+1) )
                        zcalc = rvl * TIPSLOPE * 0.5 * dist
                        b     = zt - zetac
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      ELSEIF (IPLPOS(j+1,i,k,iz).EQ.2) THEN
                        dist  = 0.5 * ( DELR(j) + DELR(j+1) )
                        zcalc = rvl * TOESLOPE * dist
                        b     = zetac - zb
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      END IF
                    END IF LCOL
C
C---------------------CHECK TIPS AND TOES FOR EACH ROW BUT THE FIRST AND LAST
                    LROW: IF ( (i.NE.1) .AND. (i.NE.NROW) ) THEN
C-----------------------BACK FACE
                      IF (IPLPOS(j,i-1,k,iz).EQ.1) THEN
                        dist  = 0.5 * ( DELC(i) + DELC(i-1) )
                        zcalc = rvl * TIPSLOPE * dist
                        b     = zt - zetac
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      ELSEIF (IPLPOS(j,i-1,k,iz).EQ.2) THEN
                        dist  = 0.5 * ( DELC(i) + DELC(i-1) )
                        zcalc = rvl * TOESLOPE * dist
                        b     = zetac - zb
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      END IF
C-----------------------FRONT FACE
                      IF (IPLPOS(j,i+1,k,iz).EQ.1) THEN
                        dist  = 0.5 * ( DELC(i) + DELC(i+1) )
                        zcalc = rvl * TIPSLOPE * dist
                        b     = zt - zetac
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      ELSE IF (IPLPOS(j,i+1,k,iz).EQ.2) THEN
                        dist  = 0.5 * ( DELC(i) + DELC(i+1) )
                        zcalc = rvl * TOESLOPE * dist
                        b     = zetac - zb
                        IF ( b.GT.zcalc ) THEN
                          IADPTMOD = 1
                          ADPTVAL  = b / zcalc
                        END IF
                      END IF
                    END IF LROW
                  END IF LACTIVE
                END DO JTIPTOE
              END DO ITIPTOE
            END DO IZTIPTOE
          END DO ZTIPTOE
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWICHKSLOPE
C
C---------MOVE TIPS AND TOES
C---------NOTE: MAY DIFFER FROM MATLAB PROTOTYPE BECAUSE ZETA IS
C---------ALTERED DURING THE SEARCH PROCESS
      SUBROUTINE SSWITIPTOE()
C-----VERSION 3 26SEP2010 SSWITIPTOE
C     ******************************************************************
C     MOVE TIPS AND TOES FOR SWI PACKAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
        USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IFREFM,
     2                        LBOTM,BOTM,
     3                        CR,CC,CV,HCOF,RHS,
     4                        DELR,DELC,IBOUND,HNEW,HOLD,
     5                        BUFF,ISSFLG,NSTP
        USE GWFBASMODULE, ONLY: DELT,HDRY,TOTIM,IHDDFL,IBUDFL
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k, n
        INTEGER :: iz, iz2, iz3
        INTEGER :: izrev
        INTEGER :: icount
        REAL :: z, zt, zb, zetac, zcalc, zdiff, zetaavg, qztop
C.........TESTING ONLY MOVING TIPS AND TOES WHEN A POSITIVE ZETAMIN IS SPECIFIED
        MOVETIPTOE: IF (ZETAMIN.GT.0.0) THEN
          ZTIPTOE: DO k = 1, NLAY
            IZTIPTOE: DO iz = 2, NZONES
              ITIPTOE: DO i = 1, NROW
                JTIPTOE: DO j = 1, NCOL
C-------------------SET TEMPORARY ZETA VARIABLES
                  zt    = ZETA(j,i,k,1)
                  zb    = ZETA(j,i,k,NZONES+1)
                  zetac = ZETA(j,i,k,iz)
                  qztop = 0.0
!                  IF ( k.GT.1 ) THEN
                  IF ( k.NE.1 ) THEN
                    qztop = CV(j,i,k-1) * (HNEW(j,i,k)-HNEW(j,i,k-1)) +
     2                      QZEXTRA(j,i,k)
                  END IF
C...................TESTLAYER
                  IF( NOPT.LT.0 )THEN
                    IF ( IPLPOS(j,i,k,iz).EQ.1 ) THEN
!                      IF (k.NE.1) THEN
                      IF (IPLPOS(j,i,k-1,iz).EQ.2) THEN
                        IF ( qztop.GE.0 ) THEN
                          ZETA(j,i,k-1,iz) =
     2                     ZETA(j,i,k-1,NZONES+1) + ZETAMIN
                        END IF
                      END IF
!                      END IF
                    END IF
                  END IF
C...................TESTLAYER
C
C-------------------ONLY EVALUATE TIP AND TOE MOVEMENT TO ADJACENT CELLS FOR CELLS WHERE
C                   ZETA SURFACE IS NOT AT THE TOP (ILPOS=1) OR BOTTOM (ILPOS=2) OF THE LAYER
                  LACTIVE: IF (IPLPOS(j,i,k,iz).EQ.0) THEN
C--------------------ADJUST TIPS AND TOES FOR EACH COLUMN BUT THE FIRST AND LAST
                    LCOL: IF ((j.NE.1).AND.(j.NE.NCOL)) THEN
C----------------------LEFT FACE
                      IF (IPLPOS(j-1,i,k,iz).EQ.1) THEN
                        zcalc = TIPSLOPE * 0.5 * (DELR(j)+DELR(j-1))
                        IF ((zt-zetac).GT.zcalc) THEN
                          IF (qztop.LE.0) THEN
                            ZETA(j,i,k,iz)=zetac +
     2                       DELZETA*DELR(j-1)/DELR(j)
                            ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,1)-DELZETA
                          END IF
                        ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
                          IF (IPLPOS(j+1,i,k,iz).EQ.0) THEN
                            ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,iz)-(zt-zetac)
                            ZETA(j,i,k,iz)=zt
                          END IF
                        END IF
                      ELSEIF (IPLPOS(j-1,i,k,iz).EQ.2) THEN
                        zcalc = zb + TOESLOPE * 0.5*(DELR(j)+DELR(j-1))
                        IF (zetac.GT.zcalc) THEN
                          ZETA(j,i,k,iz)   = zetac -
     2                                       DELZETA*DELR(j-1)/DELR(j)
                          ZETA(j-1,i,k,iz) = ZETA(j-1,i,k,NZONES+1) +
     2                                       DELZETA
                        ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
                          IF (IPLPOS(j+1,i,k,iz).EQ.0) THEN
                            ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,iz)+(zetac-zb)
                            ZETA(j,i,k,iz)=zb
                          END IF
                        END IF
                      END IF
C----------------------RIGHT FACE
                      IF (IPLPOS(j+1,i,k,iz).EQ.1) THEN
                        zcalc = TIPSLOPE * 0.5 * (DELR(j)+DELR(j+1))
                        IF ((zt-zetac).GT.zcalc) THEN
                          IF (qztop.LE.0) THEN
                            ZETA(j,i,k,iz) = zetac +
     2                       DELZETA*DELR(j+1)/DELR(j)
                            ZETA(j+1,i,k,iz)=ZETA(j+1,i,k,1)-DELZETA
                          END IF
                        ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
                          IF (IPLPOS(j-1,i,k,iz).EQ.0) THEN
                            ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,iz)-(zt-zetac)
                            ZETA(j,i,k,iz)=zt
                          END IF
                        END IF
                      ELSEIF (IPLPOS(j+1,i,k,iz).EQ.2) THEN
                        zcalc = zb + TOESLOPE * 0.5*(DELR(j)+DELR(j+1))
                        IF (zetac.GT.zcalc) THEN
                          ZETA(j,i,k,iz)   = zetac -
     2                                       DELZETA*DELR(j+1)/DELR(j)
                          ZETA(j+1,i,k,iz) = ZETA(j+1,i,k,NZONES+1) +
     2                                       DELZETA
                        ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
                          IF (IPLPOS(j-1,i,k,iz).EQ.0) THEN
                            ZETA(j-1,i,k,iz)=ZETA(j-1,i,k,iz)+(zetac-zb)
                            ZETA(j,i,k,iz)=zb
                          END IF
                        END IF
                      END IF
                    END IF LCOL
C
C---------------------ADJUST TIPS AND TOES FOR EACH ROW BUT THE FIRST AND LAST
                    LROW: IF ( (i.NE.1) .AND. (i.NE.NROW) ) THEN
C-----------------------BACK FACE
                      IF (IPLPOS(j,i-1,k,iz).EQ.1) THEN
                        zcalc = TIPSLOPE * 0.5 * (DELC(i)+DELC(i-1))
                        IF ((zt-zetac).GT.zcalc) THEN
                          ZETA(j,i,k,iz)=zetac+DELZETA*DELC(i-1)/DELC(i)
                          ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,1)-DELZETA
                        ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
                          IF (IPLPOS(j,i+1,k,iz).EQ.0) THEN
                            ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,iz)-(zt-zetac)
                            ZETA(j,i,k,iz)=zt
                          END IF
                        END IF
                      ELSEIF (IPLPOS(j,i-1,k,iz).EQ.2) THEN
                        zcalc = zb + TOESLOPE*0.5*(DELC(i)+DELC(i-1))
                        IF (zetac.GT.zcalc) THEN
                          ZETA(j,i,k,iz)   = zetac -
     2                                       DELZETA*DELC(i-1)/DELC(i)
                          ZETA(j,i-1,k,iz) = ZETA(j,i-1,k,NZONES+1) +
     2                                       DELZETA
                        ELSEIF ((zetac-zb).LT.ZETAMIN) THEN
                          IF (IPLPOS(j,i+1,k,iz).EQ.0) THEN
                            ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,iz)+(zetac-zb)
                            ZETA(j,i,k,iz)=zb
                          END IF
                        END IF
                      END IF
C-----------------------FRONT FACE
                      IF (IPLPOS(j,i+1,k,iz).EQ.1) THEN
                        zcalc = TIPSLOPE*0.5*(DELC(i)+DELC(i+1))
                        IF ((zt-zetac).GT.zcalc) THEN
                          ZETA(j,i,k,iz)=zetac+DELZETA*DELC(i+1)/DELC(i)
                          ZETA(j,i+1,k,iz)=ZETA(j,i+1,k,1)-DELZETA
                        ELSEIF ((zt-zetac).LT.ZETAMIN) THEN
                          IF (IPLPOS(j,i-1,k,iz).EQ.0) THEN
                            ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,iz)-(zt-zetac)
                            ZETA(j,i,k,iz)=zt
                          END IF
                        END IF
                      ELSE IF (IPLPOS(j,i+1,k,iz).EQ.2) THEN
                        zcalc = zb + TOESLOPE * 0.5*(DELC(i)+DELC(i+1))
                        IF (zetac.GT.zcalc) THEN
                          ZETA(j,i,k,iz)   = zetac -
     2                                       DELZETA*DELC(i+1)/DELC(i)
                          ZETA(j,i+1,k,iz) = ZETA(j,i+1,k,NZONES+1) +
     2                                       DELZETA
                        ELSE IF ((zetac-zb).LT.ZETAMIN) THEN
                          IF (IPLPOS(j,i-1,k,iz).EQ.0) THEN
                            ZETA(j,i-1,k,iz)=ZETA(j,i-1,k,iz)+(zetac-zb)
                            ZETA(j,i,k,iz)=zb
                          END IF
                        END IF
                      END IF
                    END IF LROW
                  END IF LACTIVE
                END DO JTIPTOE
              END DO ITIPTOE
            END DO IZTIPTOE
          END DO ZTIPTOE
C
C-----------CHECK WHETHER ANYWHERE THE THICKNESS GETS TOO THIN
          KZTT: DO k=1,NLAY
            IZTT: DO i=1,NROW
              JZTT: DO j=1,NCOL
                IZZTT: DO iz=2,NZONES
                  IF (IPLPOS(j,i,k,iz).EQ.0) THEN
C                    IF ((ZETA(j,i,k,iz)-ZETA(j,i,k,NZONES+1)).LT.ZETAMIN) THEN
                    zdiff = ZETA(j,i,k,iz)-ZETA(j,i,k,NZONES+1)
                    IF (zdiff.LT.ZETAMIN) THEN
                      ZETA(j,i,k,iz)=ZETA(j,i,k,NZONES+1)
                    END IF
C                    IF ((ZETA(j,i,k,1)-ZETA(j,i,k,iz)).LT.ZETAMIN) THEN
                    zdiff = ZETA(j,i,k,1)-ZETA(j,i,k,iz)
                    IF (zdiff.LT.ZETAMIN) THEN
                      ZETA(j,i,k,iz) = ZETA(j,i,k,1)
                    END IF
                  END IF
                END DO IZZTT
              END DO JZTT
            END DO IZTT
          END DO KZTT
C
C-----------MODIFY ZETA ANYWHERE THE PLANES ARE CROSSING
          KPX: DO k=1,NLAY
            IPX: DO i=1,NROW
              JPX: DO j=1,NCOL
                IZPX: DO iz=2,NZONES-1
                  IF ((ZETA(j,i,k,iz)-ZETA(j,i,k,iz+1)).LT.0.001) THEN
                    zetaavg = 0.5 * (ZETA(j,i,k,iz)+ZETA(j,i,k,iz+1))
                    ZETA(j,i,k,iz)=zetaavg
                    ZETA(j,i,k,iz+1)=zetaavg
                    IZ2PX: DO iz2=2,iz-1
                      izrev = iz+1 - iz2
                      zdiff = ZETA(j,i,k,izrev)-ZETA(j,i,k,iz+1)
                      IF (zdiff.LT.0.001) THEN
                        zetaavg = 0.
                        icount = 0
                        DO iz3 = izrev, iz+1
                          icount = icount + 1
                          zetaavg = zetaavg + ZETA(j,i,k,iz3)
                        END DO
                        IF ( icount.NE.0 ) THEN
                          zetaavg = zetaavg / icount
                          DO iz3 = izrev, iz+1
                            ZETA(j,i,k,iz3) = zetaavg
                          END DO
                        END IF
                      END IF
                    END DO IZ2PX
                  END IF
                END DO IZPX
              END DO JPX
            END DO IPX
          END DO KPX
        END IF MOVETIPTOE ! this loop done only when zetamin > 0
C
C---------RETURN
        RETURN
      END SUBROUTINE SSWITIPTOE

C-------DEALLOCATE SWI DATA FOR A GRID
      SUBROUTINE GWF2SWI1DA(Igrid)
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Igrid
C       + + + CODE + + +
        DEALLOCATE(GWFSWIDAT(Igrid)%NPLN)
        DEALLOCATE(GWFSWIDAT(Igrid)%ISTRAT)
        DEALLOCATE(GWFSWIDAT(Igrid)%NZONES)

        DEALLOCATE(GWFSWIDAT(Igrid)%NADPTFLG)
        DEALLOCATE(GWFSWIDAT(Igrid)%NADPTMX)
        DEALLOCATE(GWFSWIDAT(Igrid)%NADPTMN)
        DEALLOCATE(GWFSWIDAT(Igrid)%ADPTFCT)
        DEALLOCATE(GWFSWIDAT(Igrid)%IADPT)
        DEALLOCATE(GWFSWIDAT(Igrid)%IADPTMOD)
        DEALLOCATE(GWFSWIDAT(Igrid)%ADPTVAL)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWIDELT)

        DEALLOCATE(GWFSWIDAT(Igrid)%NOBS)
        DEALLOCATE(GWFSWIDAT(Igrid)%IOBSHEADER)

        DEALLOCATE(GWFSWIDAT(Igrid)%ISWIZT)
        DEALLOCATE(GWFSWIDAT(Igrid)%ISWICB)
        DEALLOCATE(GWFSWIDAT(Igrid)%ISWIOBS)

        DEALLOCATE(GWFSWIDAT(Igrid)%NOPT)

        DEALLOCATE(GWFSWIDAT(Igrid)%NLAYSWI)
C         SOLVER
        DEALLOCATE(GWFSWIDAT(Igrid)%NSOLVER)
!        DEALLOCATE(GWFSWIDAT(Igrid)%ZCLOSE)
        DEALLOCATE(GWFSWIDAT(Igrid)%IPRSOL)
        DEALLOCATE(GWFSWIDAT(Igrid)%MUTSOL)
C         SWI PARAMETERS
        DEALLOCATE(GWFSWIDAT(Igrid)%TOESLOPE)
        DEALLOCATE(GWFSWIDAT(Igrid)%TIPSLOPE)
        DEALLOCATE(GWFSWIDAT(Igrid)%ZETAMIN)
        DEALLOCATE(GWFSWIDAT(Igrid)%DELZETA)
        DEALLOCATE(GWFSWIDAT(Igrid)%IBO)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWIHCOF)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWISOLVCR)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWISOLVCC)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWISOLVCV)
        DEALLOCATE(GWFSWIDAT(Igrid)%ZETA)
        DEALLOCATE(GWFSWIDAT(Igrid)%ZETAOLD)
!        DEALLOCATE(GWFSWIDAT(Igrid)%ZETAIT)
        DEALLOCATE(GWFSWIDAT(Igrid)%SSZ)
        DEALLOCATE(GWFSWIDAT(Igrid)%EPS)
        DEALLOCATE(GWFSWIDAT(Igrid)%NUS)
        DEALLOCATE(GWFSWIDAT(Igrid)%DELNUS)
        DEALLOCATE(GWFSWIDAT(Igrid)%NUPLANE)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWICR)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWICC)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWICUMCR)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWICUMCC)
        DEALLOCATE(GWFSWIDAT(Igrid)%NUTOP)
        DEALLOCATE(GWFSWIDAT(Igrid)%NUBOT)
        DEALLOCATE(GWFSWIDAT(Igrid)%BRHS)
        DEALLOCATE(GWFSWIDAT(Igrid)%QZEXTRA)
        DEALLOCATE(GWFSWIDAT(Igrid)%QZEXTRACUM)
        DEALLOCATE(GWFSWIDAT(Igrid)%DUM)
        DEALLOCATE(GWFSWIDAT(Igrid)%RHSFRESH)
        DEALLOCATE(GWFSWIDAT(Igrid)%HCOFFRESH)
        DEALLOCATE(GWFSWIDAT(Igrid)%IPLPOS)
        DEALLOCATE(GWFSWIDAT(Igrid)%IZONENR)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWIOBS)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWIDE4)
        DEALLOCATE(GWFSWIDAT(Igrid)%SWIPCG)
C
C-------RETURN
        RETURN
      END SUBROUTINE GWF2SWI1DA

C
C-------SET POINTERS TO SWI DATA FOR A GRID
      SUBROUTINE SGWF2SWI1PNT(Igrid)
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Igrid
C       + + + CODE + + +
        NPLN=>GWFSWIDAT(Igrid)%NPLN
        ISTRAT=>GWFSWIDAT(Igrid)%ISTRAT
        NZONES=>GWFSWIDAT(Igrid)%NZONES

        NADPTFLG=>GWFSWIDAT(Igrid)%NADPTFLG
        NADPTMX=>GWFSWIDAT(Igrid)%NADPTMX
        NADPTMN=>GWFSWIDAT(Igrid)%NADPTMN
        ADPTFCT=>GWFSWIDAT(Igrid)%ADPTFCT
        IADPT=>GWFSWIDAT(Igrid)%IADPT
        IADPTMOD=>GWFSWIDAT(Igrid)%IADPTMOD
        ADPTVAL=>GWFSWIDAT(Igrid)%ADPTVAL
        SWIDELT=>GWFSWIDAT(Igrid)%SWIDELT

        NOBS=>GWFSWIDAT(Igrid)%NOBS
        IOBSHEADER=>GWFSWIDAT(Igrid)%IOBSHEADER

        ISWIZT=>GWFSWIDAT(Igrid)%ISWIZT
        ISWICB=>GWFSWIDAT(Igrid)%ISWICB
        ISWIOBS=>GWFSWIDAT(Igrid)%ISWIOBS

        NOPT=>GWFSWIDAT(Igrid)%NOPT

        NLAYSWI=>GWFSWIDAT(Igrid)%NLAYSWI
C         SOLVER
        NSOLVER=>GWFSWIDAT(Igrid)%NSOLVER
!        ZCLOSE=>GWFSWIDAT(Igrid)%ZCLOSE
        IPRSOL=>GWFSWIDAT(Igrid)%IPRSOL
        MUTSOL=>GWFSWIDAT(Igrid)%MUTSOL
C         SWI PARAMETERS
        TOESLOPE=>GWFSWIDAT(Igrid)%TOESLOPE
        TIPSLOPE=>GWFSWIDAT(Igrid)%TIPSLOPE
        ZETAMIN=>GWFSWIDAT(Igrid)%ZETAMIN
        DELZETA=>GWFSWIDAT(Igrid)%DELZETA
        IBO=>GWFSWIDAT(Igrid)%IBO
        SWIHCOF=>GWFSWIDAT(Igrid)%SWIHCOF
        SWISOLVCR=>GWFSWIDAT(Igrid)%SWISOLVCR
        SWISOLVCC=>GWFSWIDAT(Igrid)%SWISOLVCC
        SWISOLVCV=>GWFSWIDAT(Igrid)%SWISOLVCV
        ZETA=>GWFSWIDAT(Igrid)%ZETA
        ZETAOLD=>GWFSWIDAT(Igrid)%ZETAOLD
!        ZETAIT=>GWFSWIDAT(Igrid)%ZETAIT
        SSZ=>GWFSWIDAT(Igrid)%SSZ
        EPS=>GWFSWIDAT(Igrid)%EPS
        NUS=>GWFSWIDAT(Igrid)%NUS
        DELNUS=>GWFSWIDAT(Igrid)%DELNUS
        NUPLANE=>GWFSWIDAT(Igrid)%NUPLANE
        SWICR=>GWFSWIDAT(Igrid)%SWICR
        SWICC=>GWFSWIDAT(Igrid)%SWICC
        SWICUMCR=>GWFSWIDAT(Igrid)%SWICUMCR
        SWICUMCC=>GWFSWIDAT(Igrid)%SWICUMCC
        NUTOP=>GWFSWIDAT(Igrid)%NUTOP
        NUBOT=>GWFSWIDAT(Igrid)%NUBOT
        BRHS=>GWFSWIDAT(Igrid)%BRHS
        QZEXTRA=>GWFSWIDAT(Igrid)%QZEXTRA
        QZEXTRACUM=>GWFSWIDAT(Igrid)%QZEXTRACUM
        DUM=>GWFSWIDAT(Igrid)%DUM
        RHSFRESH=>GWFSWIDAT(Igrid)%RHSFRESH
        HCOFFRESH=>GWFSWIDAT(Igrid)%HCOFFRESH
        IPLPOS=>GWFSWIDAT(Igrid)%IPLPOS
        IZONENR=>GWFSWIDAT(Igrid)%IZONENR
        SWIOBS=>GWFSWIDAT(Igrid)%SWIOBS
        SWIDE4=>GWFSWIDAT(Igrid)%SWIDE4
        SWIPCG=>GWFSWIDAT(Igrid)%SWIPCG
C
C---------RETURN
        RETURN
      END SUBROUTINE SGWF2SWI1PNT
C
C-------SAVE POINTERS TO SWI DATA FOR A GRID
      SUBROUTINE SGWF2SWI1PSV(Igrid)
        USE GWFSWIMODULE
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: Igrid
C       + + + CODE + + +
        GWFSWIDAT(Igrid)%NPLN=>NPLN
        GWFSWIDAT(Igrid)%ISTRAT=>ISTRAT
        GWFSWIDAT(Igrid)%NZONES=>NZONES

        GWFSWIDAT(Igrid)%NADPTFLG=>NADPTFLG
        GWFSWIDAT(Igrid)%NADPTMX=>NADPTMX
        GWFSWIDAT(Igrid)%NADPTMN=>NADPTMN
        GWFSWIDAT(Igrid)%ADPTFCT=>ADPTFCT
        GWFSWIDAT(Igrid)%IADPT=>IADPT
        GWFSWIDAT(Igrid)%IADPTMOD=>IADPTMOD
        GWFSWIDAT(Igrid)%ADPTVAL=>ADPTVAL
        GWFSWIDAT(Igrid)%SWIDELT=>SWIDELT

        GWFSWIDAT(Igrid)%NOBS=>NOBS
        GWFSWIDAT(Igrid)%IOBSHEADER=>IOBSHEADER

        GWFSWIDAT(Igrid)%ISWIZT=>ISWIZT
        GWFSWIDAT(Igrid)%ISWICB=>ISWICB
        GWFSWIDAT(Igrid)%ISWIOBS=>ISWIOBS

        GWFSWIDAT(Igrid)%NOPT=>NOPT

        GWFSWIDAT(Igrid)%NLAYSWI=>NLAYSWI
C         SOLVER
        GWFSWIDAT(Igrid)%NSOLVER=>NSOLVER
!        GWFSWIDAT(Igrid)%ZCLOSE=>ZCLOSE
        GWFSWIDAT(Igrid)%IPRSOL=>IPRSOL
        GWFSWIDAT(Igrid)%MUTSOL=>MUTSOL
C         SWI PARAMETERS
        GWFSWIDAT(Igrid)%TOESLOPE=>TOESLOPE
        GWFSWIDAT(Igrid)%TIPSLOPE=>TIPSLOPE
        GWFSWIDAT(Igrid)%ZETAMIN=>ZETAMIN
        GWFSWIDAT(Igrid)%DELZETA=>DELZETA
        GWFSWIDAT(Igrid)%IBO=>IBO
        GWFSWIDAT(Igrid)%SWIHCOF=>SWIHCOF
        GWFSWIDAT(Igrid)%SWISOLVCR=>SWISOLVCR
        GWFSWIDAT(Igrid)%SWISOLVCC=>SWISOLVCC
        GWFSWIDAT(Igrid)%SWISOLVCV=>SWISOLVCV
        GWFSWIDAT(Igrid)%ZETA=>ZETA
        GWFSWIDAT(Igrid)%ZETAOLD=>ZETAOLD
!        GWFSWIDAT(Igrid)%ZETAIT=>ZETAIT
        GWFSWIDAT(Igrid)%SSZ=>SSZ
        GWFSWIDAT(Igrid)%EPS=>EPS
        GWFSWIDAT(Igrid)%NUS=>NUS
        GWFSWIDAT(Igrid)%DELNUS=>DELNUS
        GWFSWIDAT(Igrid)%NUPLANE=>NUPLANE
        GWFSWIDAT(Igrid)%SWICR=>SWICR
        GWFSWIDAT(Igrid)%SWICC=>SWICC
        GWFSWIDAT(Igrid)%SWICUMCR=>SWICUMCR
        GWFSWIDAT(Igrid)%SWICUMCC=>SWICUMCC
        GWFSWIDAT(Igrid)%NUTOP=>NUTOP
        GWFSWIDAT(Igrid)%NUBOT=>NUBOT
        GWFSWIDAT(Igrid)%BRHS=>BRHS
        GWFSWIDAT(Igrid)%QZEXTRA=>QZEXTRA
        GWFSWIDAT(Igrid)%QZEXTRACUM=>QZEXTRACUM
        GWFSWIDAT(Igrid)%DUM=>DUM
        GWFSWIDAT(Igrid)%RHSFRESH=>RHSFRESH
        GWFSWIDAT(Igrid)%HCOFFRESH=>HCOFFRESH
        GWFSWIDAT(Igrid)%IPLPOS=>IPLPOS
        GWFSWIDAT(Igrid)%IZONENR=>IZONENR
        GWFSWIDAT(Igrid)%SWIOBS=>SWIOBS
        GWFSWIDAT(Igrid)%SWIDE4=>SWIDE4
        GWFSWIDAT(Igrid)%SWIPCG=>SWIPCG
C
C---------RETURN
        RETURN
      END SUBROUTINE SGWF2SWI1PSV
